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The dynamics of complex systems in nature often occurs in terms of punctuations, or avalanches, 
rather than following a smooth, gradual path. A comprehensive theory of avalanche dynamics in 
models of growth, interface depinning, and evolution is presented. Specifically, we include the Bak- 
Sneppen evolution model, the Sneppen interface depinning model, the Zaitsev flux creep model, 
invasion percolation, and several other depinning models into a unified treatment encompassing a 
large class of far from equilibrium processes. The formation of fractal structures, the appearance 
of 1// noise, diffusion with anomalous Hurst exponents, Levy flights, and punctuated equilibria 
can all be related to the same underlying avalanche dynamics. This dynamics can be represented 
as a fractal in d spatial plus one temporal dimension. The complex state can be reached either 
by tuning a parameter, or it can be self-organized (SOC). We present two exact equations for the 
avalanche behavior in the latter case. (1) The slow approach to the critical attractor, i.e. the 
process of self-organization, is governed by a "gap" equation for the divergence of avalanche sizes. 
(2) The hierarchical structure of avalanches is described by an equation for the average number of 
sites covered by an avalanche. The exponent 7 governing the approach to the critical state appears 
as a constant rather than as a critical exponent. In addition, the conservation of activity in the 
stationary state manifests itself through the superuniversal result rj = 0. The exponent n for the 
Levy flight jumps between subsequent active sites can be related to other critical exponents through 
a study of "backward avalanches." We develop a scaling theory that relates many of the critical 
exponents in this broad category of extremal models, representing different universality classes, to 
two basic exponents characterizing the fractal attractor. The exact equations and the derived set 
of scaling relations are consistent with numerical simulations of the above mentioned models. 



I. INTRODUCTION 

The term spatio-temporal complexity describes sys- 
tems that contain information over a wide range of length 
and time scales |jj . Specifically, if such a system is stud- 
ied through a magnifying glass, then no matter what the 
level of magnification is, a variation of the image is seen 
as different parts of the system are viewed. Similarly, 
if the image of a local part of the system is averaged 
over a time window, different images are seen at different 
times, no matter how large the time window is. This con- 
trasts with the behavior of ordered systems, or random 
systems, including chaotic ones, which become uniform 
when viewed at sufficiently large length or time scales. 

The appearance of complexity in nature presents a fas- 
cinating, long-standing puzzle. In this article, a qualita- 
tive and quantitative theory for the dynamics of complex 
systems is presented in the context of simple mathemat- 
ical models for biological evolution and growth phenom- 
ena far from equilibrium. Spatio-temporal complexity 
emerges as the result of avalanche dynamics in driven 
systems. We unify the origin of fractals, 1// noise, 



Hurst exponents for anomalous diffusion, Levy flights, 
and punctuated equilibrium behavior, and elucidate their 
relationships through analytical and numerical studies of 
simple models, defined in Section II. These phenomena 
are signatures of spatio-temporal complexity and are re- 
lated via scaling relations to the fractal properties of the 
avalanches, as summarized in Table |. 

We wish to focus on three general, empirical phenom- 
ena that are manifestations of complexity. First, Man- 
delbrot ^ has pointed out the widespread occurrence 
of self-similar, fractal behavior in nature. Mountains, 
coastlines, and perhaps even the Universe jU have fea- 
tures at all scales. River networks consist of streams of 
all sizes ||] ; the pattern of earthquake faults, cracks in 
rocks, or oil reservoirs is self-similar Second, noise 
with a l/f d power spectrum is emitted from a variety 
of sources, including light from quasars Q, river flow 
, and brain activity Third, many natural and so- 
cial phenomena, including earthquakes, economic activ- 
ity, and biological evolution appear to evolve intermit- 
tently, with bursts, or avalanches extending over a wide 
range of magnitudes, rather than smoothly and gradu- 
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ally. For instance, the distribution of earthquake magni- 
tudes obeys the Gutenberg- Richter law ||] . Recent exper- 
iments on slowly driven sandpiles and rice piles 
show avalanches, or disturbances in the heap, of all sizes. 
Field et al monitored the dynamics of superconducting 
vortices and found the superconducting analog of sand- 
pile avalanches to be a power law over two decades [[l2| . 
Gould and Eldredge have proposed that biological evo- 
lution takes place in terms of punctuations, where many 
species become extinct and new species emerge, inter- 
rupting quiet periods of apparent equilibrium, known as 
stasis p3[ . Large events, such as massive extinctions or 
even the large scale structure of the Universe, may have a 
statistical explanation as the power-law tail of a Levy dis- 
tribution. This is the case, even though the large events 
are identifiable and tend to be viewed, erroneously we 
believe, as "atypical". 

Even though spatio-temporal complexity is ubiquitous 
in nature, until recently, little understanding of its ori- 
gin has been achieved. One clear exception, though, are 
critical points for second order phase transitions in equi- 
librium systems, which are characterized by scale invari- 
ance. Borrowing the terminology from equilibrium ther- 
modynamics, we shall also refer to systems with a large 
range of length and time scales as "critical." Spatial 
complexity occurs at the percolation transition in ran- 
dom bond or site models O] ; temporal complexity exists 
at the transition to chaos in the Feigenbaum map [ p5[ . 
These systems are critical, and thus complex, as the re- 
sult of tuning: the temperature, or some other parame- 
ter, is set to a specific value in order to achieve criticality. 
In nature, though, fine tuning of specific parameters is 
rare and unlikely; in addition, many systems in nature 
are far from an equilibrium state. 

Thus another mechanism has been proposed; namely 
systems that are far from equilibrium become critical 
through self-organization fl6| . They evolve through tran- 
sient states, which are not critical, to a dynamical attrac- 
tor poised at criticality. In order for self-organization to 
occur, these systems must be driven slowly through a 
succession of metastable states. The system jumps from 
one metastable state to another by avalanche dynamics. 
These avalanches build up long range correlations in the 
system. Here, we shall be mainly concerned with the self- 
organized critical (SOC) origin of spatio-temporal com- 
plexity. In some cases, though, similar considerations, 
such as scaling relations, can also be applied in cases 
where the criticality is tuned rather than self-organized. 
In particular, this is relevant for systems which undergo 
dcpinning transitions when a parameter is varied, such 
as interfaces, charge density wave systems, and super- 
conducting flux lattices. In this case, too, long lived 
metastable states are important and the dynamics is that 
of avalanches p 7 ?] . Avalanche dynamics in these tuned 
depinning transitions are discussed in Section VI. 

SOC complements the concept of chaos, where simple 



systems with few degrees of freedom can display com- 
plicated behavior. Chaos is associated with a fractal 
"strange" attractor in phase space. These self-similar 
structures have little to do with fractals in large spa- 
tially extended systems that have many degrees of free- 
dom. In addition, chaotic systems exhibit white noise 
(e.g. a positive Lyaponov exponent) with limited tempo- 
ral correlations whereas complex systems have long range 
spatio-temporal correlations. 

From the earliest BTW sandpile models [|l6| and, later, 
earthquake models |l8|], most of the evidence for SOC be- 
havior has been numerical. Exceptions include the work 
on singular diffusion by Carlson, Chayes, Grannan, and 
Swindle | ft9f , and the one dimensional forest fire model 
|f20f , where exact results have been found by Drossel, 
Clar, and Schwabl and also by Paczuski and Bak 
|p2| . Dhar (2^] was able to characterize the critical attrac- 
tor of the BTW sandpile model in terms of an Abelian 
group, and thereby calculate the number of states on 
the attractor. With Ramaswamy, he solved a directed 
sandpile model exactly |24|]. So far, though, none of the 
solutions of these different models have yielded a general 
phenomenology for SOC behavior. Most importantly, the 
fundamental mechanism for the self-organization process, 
via avalanches, has not been well understood. Pietronero 
and collaborators have developed a scheme, the Fixed 
Scale Transformation ^5|, and applied it to a variety 
of nonequilibrium dynamical systems including Diffusion 
Limited Aggregation p6j, the sandpile models ^7j, and 
the Bak-Sneppen evolution model p8[ . Here, we take 
a different approach where, based on certain exact re- 
sults together with a scaling ansatz, we develop a theory 
that relates different critical exponents, including the ap- 
proach to the attractor, to two basic exponents which are 
model dependent. 

A common feature of many models exhibiting SOC is 
the selection of extremal sites for the initiation of events, 
rather than statistically typical sites. This feature of ex- 
tremal dynamics has been somewhat obscured by the in- 
troduction of models, such as the BTW sandpile model 
or the forest fire model, which appear to be driven ran- 
domly. In the "random" BTW model, sand is added to 
random sites until a local threshold is exceeded and a 
toppling occurs. However, the statistics of avalanches in 
the BTW model is not changed in "deterministic" ver- 
sions where the height of every site is raised uniformly 
until one site, the least stable site, topples. In the latter 
case, randomness only enters into the initial conditions. 
One might say that in the "random" BTW model, the 
system selects the extremal sites itself, while in the "de- 
terministic" case, it is forced to do so. Similarly, in the 
Olami, Feder, Christensen earthquake model Jig] , the 
force is raised uniformly until the site with the largest 
force ruptures. Zaitsev p9| ] introduced a model for low 
temperature flux creep where the motion always takes 
place at the site with the largest force, and showed that 
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the model self-organizes to a critical state. In the forest 
fire model, it can be shown rigorously that the process is 
driven by the burning of the largest forests with the old- 
est trees, despite the fact that the trees grow randomly 
pl| , p2| . Again, burning trees randomly only serves as a 
vehicle for the system to burn the largests forests. In 
a remarkable paper, Miller et al |30) suggested through 
a different line of thinking that 1// noise and fractal- 
ity are related to the preferential selection and thermal 
activation of extremal sites. 

Recently, a variety of simple models with extremal dy- 
namics that exhibit SOC have been introduced. These 
models, including the Bak-Sneppen evolution model 
the Sneppen interface depinning model |3^ | , and the Zait- 
sev model [^9|, are defined in Section II. They are driven 
by sequentially updating the site with globally extremal 
values of the signal. This information is propagated 
through the system via local interactions. These mod- 
els, representing different universality classes, are nev- 
ertheless similar to invasion percolation |3^,Q. Inter- 
estingly, invasion percolation was recognized as an SOC 
phenomenon by Roux and Guyon p5| some years ago. 
These authors, inspired by the sandpile analogy , de- 
fined avalanches in invasion percolation in terms of a 
fluctuating extremal signal. Some additional new results 
for invasion percolation have recently been presented by 
Maslov ]|(| using methods similar to those applied in this 
work. 

Our main analytical results for this class of models, 
including invasion percolation, are encapsulated by two 
exact equations plus a superuniversal scaling law for sta- 
tionary critical processes. These results address basic 
questions that arise in SOC. Perhaps the most obvious 
question is how the system self-organizes into the com- 
plex, critical state. This is described by a "gap" equation 
that relates the rate of approach to the attractor to the 
average avalanche size |37|. This equation demonstrates 
that the stationary state of the system is a critical state 
for the avalanches, where the average avalanche size di- 
verges. Assuming that this divergence has a character- 
istic exponent 7, we show from the gap equation that 
the system approaches its attractor algebraically with a 
characteristic exponent, p = 1/(7 — 1), for the transient. 

The off-critical exponent for the transient, 7, can be 
expressed in terms of the exponents in the stationary 
state itself, via a "gamma" equation for the hierarchy of 
avalanches |||. This may result from the fact that the 
off-critical direction is stable (in SOC the critical point 
is an attractor for the dynamics) rather than being un- 
stable. Lastly, in SOC the critical point is constrained 
to be stationary. This leads to a fundamental law for 
the critical states; 77 = in all dimensions pjj. One 
dramatic consequence of this law is that the fractal di- 
mension of the active sites, df, is fixed by the probability 
distribution for avalanche sizes in the stationary state, 
i.e. df = D(t — 1). Thus the widespread appearance 



of fractal structures can, perhaps, be viewed as a conse- 
quence of the existence of a stationary limit. 

By studying the space-time behavior of the activity 
pattern in the critical state, i.e. the spatial location of 
the extremal site at a particular point in time, one can 
describe the activity pattern as a fractal embedded in 
d spatial dimensions plus one temporal dimension fl40| . 
This fractal has a mass dimension, or avalanche dimen- 
sion D. Long range time correlations, e.g. 1// noise, and 
spatial fractal behavior are unified as different cuts in this 
underlying space-time fractal. The temporal evolution at 
a specific position is expressed as the activity along a one 
dimensional time line piercing the fractal perpendicular 
to the spatial dimensions. The fractal spatial structure is 
found by cutting the fractal along the spatial directions 
at that time. We establish a formal relation between 1/ / 
type noise and fractal spatial behavior in terms of the 
avalanche dimension D. In the critical state, the dynam- 
ics is best characterized in terms of scale-free avalanches, 
initiated at extremal sites, and propagating by an anoma- 
lous diffusion process. Fig. [I] shows this space time frac- 
tal for the one dimensional evolution model. 

We have studied, in more detail, the value of the ex- 
tremal signal in time. Time directed avalanches are nat- 
urally defined in terms of this fluctuating signal |S6| . 
These avalanches have a hierarchical structure of val- 
leys within valleys. Forward and backward time directed 
avalanches have different statistical distributions in the 
stationary state, reflecting the irreversibility of extremal 
processes. The distribution of all forward avalanches, 
starting at each update step for the extremal dynamics, 
is a power law with a superuniversal exponent, rjr = 2. 
The distribution of all backward avalanches is also a 
power law, but with a different model dependent expo- 
nent, t£ u = 3 — t, where r is the usual avalanche size 
distribution exponent. 

Taken together, these considerations lead to many scal- 
ing relations for various physical quantities. All of the 
critical exponents that we consider can be expressed in 
terms of two basic exponents, for instance, the avalanche 
dimension D, and r, which characterizes the distribution 
of avalanche sizes. These scaling relations are summa- 
rized in Table [| They provide numerous points to test 
theoretical predictions with numerical simulations of dif- 
ferent models. We have made numerical tests of essen- 
tially all the scaling relations for many of the models and 
find a pattern of consistency which confirms the predicted 
scaling relations across different universality classes; nev- 
ertheless more accurate tests are needed for any specific 
result. The results of our simulations are presented in 
Table g and the error bars represent statistical errors 
for system sizes studied. We urge that extensive, accu- 
rate simulations be performed. Indeed, others |^], fL^ ] 
have already provided further confirmation of our scaling 
relations. 

It is important to note that for some models, the crit- 
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ical point can be reached in different ways. This is espe- 
cially clear in the context of depinning. The depinning 
transition can be self-organized or reached by tuning ei- 
ther an external driving force or the velocity. In the 
context of interface depinning and invasion percolation, 
the SOC version corresponds to driving at constant ve- 
locity in the limit where the velocity vanishes. Some of 
the exponents are the same in the different cases but oth- 
ers, in general, are different. For example, 77 = for all 
SOC depinning models, but may be non-zero for the de- 
pinning transition at constant driving force. The critical 
points that are reached by the self-organizing process are 
different than the critical points that are reached in an 
equivalent model by tuning a driving force. Thus, de- 
spite the fact that these differences disappear in mean 
field theory, SOC cannot, even in principle, be viewed as 
sweeping a system through a critical point, in contrast to 
the claims in Ref. [ft3[ . The similarities and differences 
between constant force and SOC depinning are elabo- 
rated in Section VI. Our scaling relations are compared 
with recent numerical simulations for constant force de- 
pinning @, @. 

Self-organized fractal growth is fundamentally different 
than growth processes that are, for example, described 
by (variants of) the Kardar-Parisi-Zhang (KPZ) equation 
p6| , The KPZ equation is scale invariant by sym- 

metry, and thus the criticality is not self-organized. It is 
essential to have long lived, metastable states for the self- 
organization process to take place through avalanches, 
without reverting to a "ground state" at or near equi- 
librium. In addition, unlike "generic scale invariance" 
p7j |, SOC does not require local conservation laws. With 
the exception of the Zaitsev model, all of the models 
we consider here are nonconservative; in spite of this, 
they can be shown to self-organize to a critical state. 
In the context of interface growth, the dynamic scaling 
approach has been used to separate avalanche dynamics 
associated with SOC and Langevin dynamics associated 
generic scale invariance into distinct phenomenological 
categories Q. 

Here, we present a comprehensive and detailed account 
of our work on the Sneppen interface model, the Bak- 
Sneppen evolution model, invasion percolation, and other 
SOC (and non SOC) critical models, including interface 
depinning. Some of these results have been previously 
published in short accounts. For clarity, we provide here 
a complete, self-contained discussion of these models and 
our most accurate and extensive numerical results. 

Section II introduces the models for the general reader. 
Section III examines the transient self-organization pro- 
cess and introduces the "gap" equation. Section IV dis- 
cusses the avalanche hierarchy in the stationary state. In 
particular we present the "gamma" equation, the law for 
stationary states, r\ = 0, and a discussion of time-directed 
avalanches. The concept of backward avalanches is used 
to determine the exponent 7r governing the distribution 



of jumps in the activity, which has a Levy distribution. 
Section V unifies spatial fractal behavior with 1// type 
noise. Section VI contains our results on interface depin- 
ning. In the concluding section, the new scaling relations 
that we derive are summarized in Table B, and our nu- 
merical results are summarized in Table ffl . Appendix A 
explains in more detail the results for invasion percola- 
tion. 



II. DEFINITION OF THE MODELS 



In this section, we define all of the models studied here. 



Evolution [31]: The Bak-Sneppen evolution model is 
perhaps the simplest model of SOC. In this 'toy' 
evolution model, random numbers /; are assigned 
independently to each site on a d-dimensional lat- 
tice of linear size L. They are chosen from a uni- 
form distribution between zero and one, V(f). At 
each update step, the extremal site, i.e. the one 
with the smallest random number, f m im is located. 
That site, and its 2d nearest neighbors are then as- 
signed new random numbers, drawn independently 
from the flat distribution V . After many updates 
have occurred, the system reaches a statistically 
stationary state in which the density of random 
numbers in the system vanishes for / < f c and is 
uniform above f c . In the thermodynamic L — » 00 
limit, no random number with / > f c is ever the 
extremal site. A snapshot of the stationary state in 
a finite one dimensional system is shown in Fig. |^. 
Except for a localized region, the avalanche, where 
there are small random numbers, all the random 
numbers in the system have values above the self- 
organized threshold f c = 0.66702 ± 0.00003 in one 
dimension. 

The evolution model exhibits punctuated equilib- 
rium behavior as does real biology It simu- 
lates evolutionary activity in terms of mutations of 
individual species and interdependencies in a food 
chain. The random numbers represent the fitness of 
a species, and chosing the smallest random number 
mimicks the "Darwinian" principle that the least 
fit species is replaced or mutates. The dynamical 
impact of the mutation of the least fit species on 
the rest of the ecology is simulated by changing the 
fitness of neighboring species on the lattice. Dis- 
cussions of its possible connection to biological evo- 
lution and macroeconomics may be found in Refs. 
fl) , f£|] , and a mean field analysis in @ , ^ , [f§ . 
A generalized M -vector model where each species 
has many (M) internal degrees of freedom has very 
recently been introduced and solved exactly p3[ | in 
the M — > 00 limit. 
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Sneppen [32]: In the Sneppen interface model for SOC 
depinning, an interface of size L d defined on a dis- 
crete lattice (x, h) moves under the influence of 
quenched random pinning forces f(x, h) assigned 
independently from the flat distribution V . Ini- 
tially, h(x) — 0. Growth occurs by advancing the 
extremal site on the interface with the smallest ran- 
dom pinning force, f m i n , by one step. Then a con- 
straint is imposed for all nearest neighbor gradi- 
ents, \h(x) — h(x')\ < 1. This is met by advanc- 
ing the heights of neighboring sites. The process 
is repeated indefinitely. Like the evolution model, 
the Sneppen model also reaches a statistically sta- 
tionary state where the density of random pinning 
forces on the interface vanishes for / < f c and is 
uniform above f c . In this state, the interface moves 
in bursts of localized activity, as indicated schemat- 
ically in Fig. H Outside of these regions of activity, 
the interface is frozen for potentially long periods 
of time until a burst moves into the frozen region. 
The interface moves in a jerky, irregular manner, 
rather than smoothly and gradually advancing as a 
whole. 



represents a discretized Laplacian. Thus the local 
internal force is 



F int (x,t) = V 2 h(x,t) + f(x,h(x,t)) 



(1) 



Tang and Leschhorn 54 showed that the one- 



dimensional model, in the stationary state, iden- 
tifies from time to time with paths on a critical 
directed percolation cluster. It has been proposed 
that this identification with blocking paths can be 
used to obtain all of the critical exponents in terms 
of two independent exponents ]55[ ] . The identifica- 
tion with "blocking surfaces" is analogous to inva- 
sion percolation, where the invading region iden- 
tifies from time to time with a critical cluster of 
ordinary percolation J|3|] , |35| . Possible physical re- 
alizations of the Sneppen model and other interface 
depinning models are discussed in Section VI. 



Zaitsev [ 29 1 : Zaitsev introduced an extremal model for 
flux creep. Random numbers /j, chosen from V , are 
assigned independently to each site on a d dimen- 
sional lattice of linear size L. At each update step, 
the site with the largest number is chosen, and a 
random number chosen from V is subtracted from 
that site and equally distributed to the 2d nearest 
neighbors. The activity conserves the sum of all 
the random numbers in the system. It has been 
suggested Q that this model is in the same uni- 
versality class as a self-organized "linear" interface 
depinning model (LIM), sometimes referred to as 
the quenched Edwards- Wilkins J57|] equation. 

LIM: In the linear interface depinning model, the force 
at each site has a random contribution, f(x,h), 
chosen from V, to represent quenched random 
pinning forces, plus a linear configurational term 
fconf ~ V 2 /i, where h is the local height and V 2 



We consider cases where this model is driven either 

(a) with extremal dynamics |4C|j5^ ] or (b) in a paral- 
lel non self-organized mode [|58| . In (a) the site with 
the maximum total force is advanced by one unit, 
leading to SOC, or constant velocity depinning. In 

(b) the model may be tuned to a depinning tran- 
sition by adding an external driving force F to all 
sites and advancing in parallel every site, where the 
total force {Fi nt +F) is positive, by one step. When 
F > F c the interface moves with a finite velocity, 
while for F < F c it becomes stuck in a metastablc 
state. At F = F c it undergoes a depinning tran- 
sition. A tuned depinning transition may also be 
realized by externally setting the velocity, v, of the 
interface and allowing the force F to fluctuate so as 
to maintain that velocity. The stationary state of 
the SOC version corresponds to the depinning tran- 
sition at constant velocity. The SOC version tunes 
itself to a constant velocity depinning transition. In 
Section VI, a comparison is made between the LIM 
and models for fluid invasion in porous media J59| , 
jBOf and interface depinning in the Random Field 
Ising Model |6^]. There, it is also argued that r 
and D are the same for the (constant force) tuned 
and SOC variants, but rj and other dynamical ex- 
ponents (e.g. df and z) can be different. 

TLB: A model for depinning of an elastic interface at 
constant force that was studied by Tang-Leschhorn 
|62| and also by Buldyrev et al Again, an in- 
terface of size L d defined on a discrete lattice (x, h) 
moves under the influence of random pinning forces 
f(x, h) assigned independently from V . Initially, 
h(x) = 0. Instead of advancing the most unsta- 
ble site, as in the Sneppen model, all unstable sites 
with f < F are advanced in parallel. Then the 
constraint is imposed for all nearest neighbor gra- 
dients, \h(x) — h(x')\ < 1. The system is relaxed 
to meet the gradient constraint. When F = F c the 
interface undergoes a depinning transition. Tang 
and Leschhorn H] , 1 62 showed that in the critical 
state both the TLB model and the Sneppen model, 
in one dimension, identify with a directed percola- 
tion cluster of sites with / > / c . In Section VI, it 
is argued that the exponents r and D are the same 
for the TLB and Sneppen models; while rj, z, and 
other dynamical exponents can be different. 

DP: In directed percolation, a preferred direction, la- 
beled by t, is chosen and bonds are oriented with 
respect to t. Percolation is only allowed in the di- 
rection of increasing t. Each bond exists with prob- 
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ability /. When / = f c , the DP cluster can become 
infinitely large. This model can be viewed as the 
"parallel" constant force version of the evolution 
model. In DP, there is a critical point which is 
tuned. 



Invasion Percolation [33], [34] 



Invasion percolation is the oldest member of this 
class of models. It was studied as an SOC phe- 
nomenon, in terms of avalanches, by Roux and 
Guyon p5| |. In invasion percolation, random num- 
bers fi are assigned independently to each site on a 
d + 1-dimensional lattice of linear size L. They are 
chosen from a uniform distribution between zero 
and one, V(f). Initially, one d-dimensional side of 
the lattice is the "invaded cluster" . The random 
numbers at the boundary of the invaded cluster 
are examined. At each update step, the extremal 
site, i.e. the site with the smallest random number, 
fmin , on the boundary of the invaded cluster is lo- 
cated. That site is then invaded, and new sites sus- 
ceptible to growth are added to the boundary. The 
universality class of the model is sensitive to the ex- 
act definition of the boundary. Different definitions 
of the boundary (IP-1, IP-2, IP-3) are given in Ap- 
pendix A, along with a discussion of their physical 
realizations. 



III. SELF-ORGANIZATION 

Self-organization, as used here, refers to a dynamical 
process whereby a system starts in a state with uncor- 
related behavior and ends up in a complex state with 
a high degree of correlation. The time that it takes to 
self-organize grows as the system size increases; so that 
for large systems self-organization is a slow process. In 
contrast to the earlier models of SOC, the process of self- 
organization in the extremal models that are discussed 
here is very well understood. The stationary state is crit- 
ical and is approached algebraically, through a sequence 
of transient states. In order to demonstrate this, we shall 
first discuss self-organization in the context of the evo- 
lution model, which has been defined in Section II, and 
then generalize the results to other models. 

Let us consider the situation where the distribution 
of /'s initially is uniform in the interval [0, 1] in a d- 
dimensional evolution model of linear size L. Initially, the 
activity, i.e. the spatial location of the smallest random 
number in time, jumps randomly throughout the system. 
Eventually, after a long transient, the system reaches a 
complex state where the activity is correlated, as shown 
in Fig. ^. In order to study this transient process, we 
examine the value of the minimal random number chosen, 
fmin, as a function of sequential time, s, or the total 



number of updates. This signal f m in(s) can be related 
to the distribution of random numbers in the system. 

The first value, / m m(0), to be chosen for updating is 
0(L~ d ). Since, by definition, there are no random num- 
bers in the system smaller than / m i n (0), this quantity 
is defined to be the initial value of the gap, G, in the 
distribution of /'s; that is G(0) = f m m(0). Eventually, 
after s updates, a larger gap G(s) > G(0) opens up. 
The density of sites with random numbers below G van- 
ishes in the thermodynamic L — > oo limit, and is uniform 
above G. The current gap, G(s), is the maximum of all 
the minimum random numbers chosen, f m i n (s'), for all 
< s' < s. Fig. H shows f min as a function of s during 
the transient for a small system. The solid line shows 
the gap, G(s), as a step-wise increasing function of s. 
The gap is an envelope function that tracks the increas- 
ing peaks in f m i n . Clearly, when the gap jumps to a new 
higher value, there are no sites in the system with random 
numbers less than the gap. Since the distribution, P(f), 
that each of the random numbers are chosen from is flat, 
all of the random numbers in the system are uniformly 
distributed in the interval [G(s), 1] at those moments in 
time when the gap jumps. Thus, the envelope function 
tracks the distribution of random numbers in the sys- 
tem. At the moments in time when the gap jumps, each 
of the random numbers in the system is independently 
distributed in the interval [G(s), 1]. 

By definition, the separate instances when the gap 
G(s) jumps to its next higher value are separated by 
avalanches. Avalanches correspond to plateaus in G dur- 
ing which f m in(s) < G(s). This guarantees that the 
events within a single avalanche are causally and spa- 
tially connected. A new avalanche is initiated every time 
the gap jumps, and all the consecutive random numbers 
which are smaller than the gap after this event must have 
originated from the site that started the new avalanche. 
Once an avalanche is over it does not affect the behavior 
of any subsequent avalanche, except in terms of the gap 
threshold. As the gap increases, the probability for the 
new random numbers to fall below the gap increases also, 
and longer and longer avalanches typically occur. 

Since the distribution of random numbers above the 
gap is flat, the average size of the jump in the gap at 
the completion of each avalanche is (1 — G(s))/L d . The 
other quantity that is needed in order to describe the self- 
organizing system is the average size of the plateau for a 
given value of G, or the average avalanche size < S >g(s) • 
We shall prove below that in the limit of large system 
sizes L, the growth of the gap versus time, s, obeys the 
exact gap equation [p7|: 



dG(s) 
ds 



1 - G(s) 



L d <S> 



G(s) 



(2) 



As the gap increases, so docs the average avalanche 
size < S >g(s), which eventually diverges as G(s) — » f c . 
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Whereupon, the model is critical and the process achieves 
stationarity. In the limit L — > oo, the density of sites with 
/ < fc vanishes, and the distribution of /'s is uniform 
above f c . The gap equation (||) defines the mechanism of 
approach to the self-organized critical attractor. It con- 
tains the essential physics of SOC phenomena. When the 
average avalanche size diverges, < S >g{s)~ * °°> the sys- 
tem becomes critical. At the same time, 4jj approaches 
zero, which means that the system becomes stationary. 
Examining Eq. (|^) one notices another mathematical 
possibility for the stationary limit, which does not re- 
quire diverging avalanches. The time derivative is 
exactly zero if the current gap G(s) — 1. This situation, 
however, is not realised in any extremal model in the 
limit that L — > oo, since the presence of any interaction 
between sites (such as replacing nearest neighbors with 
new random numbers in the evolution model) will result 
in f c < 1. 

We prove the gap equation as follows: for any selected 
resolution AG <C 1 along the gap axis there is a system 
size L sufficiently large that many avalanches are needed 
to increase the gap from G(s) to G(s) + AG. In this 
case, the sum of (a) the temporal durations S of an in- 
dividual avalanche and (b) the jumps in the gap at the 
end of each avalanche will both average within this inter- 
val because they are both sums of independent random 
variables. Therefore, in the limit L — > oo Eq. (^) is ex- 
act. To be more specific, suppose that the current value 
of the gap in the system is G(s). The average num- 
ber of avalanches required to increase the gap by AG 
is N av — AGL d /(l — G(s)). By selecting system size 
L 3> AG _1 / d we ensure that N av ^> 1. N av can be made 
arbitrarily large by taking the large L limit. In this limit, 
the average number of time steps required to increase 
the gap from G(s) to G(s) + AG is given by the interval 
As = {S) G(s) N av = {S) G(s) AGL d /{\ - G(s)), and from 
the law of large numbers the fluctuations of this inter- 
val around its average value vanish. This equation shows 
that the ratio of the interval length As to the time s re- 
quired to reach G(s), As/s, vanishes as AG — > 0. Rewrit- 
ing the interval equation as = — an d taking 
the continuum limit we restore the differential equation 




In order to solve the gap equation we need to know pre- 
cisely how the average avalanche size < S >g(s) diverges 
as the critical state is approached. It is natural to assume 
that this divergence has a power law form. In analogy 
with percolation clusters, we now make an ansatz that 
the average avalanche size diverges as G approaches f c 
as 

<S>~(f c -G)-< . (3) 

We will show in the next section that for the extremal 
models studied here the exponent 7 > 1. Let us first 
consider the case 7 > 1. Inserting Eq. (||) into Eq. (||) 



and integrating, one obtains 
A/ = / c -G(s)~(-^)-", where p = -i-j . (4) 

Eq. (Q) was also found by Ray and Jan It shows 

that the critical point (A/ = 0) is approached alge- 
braically with an exponent — — Eq. ([5]) holds over 

the range L d <C s <C L D . The lower limit requires 
that the avalanches are in the scaling regime, so that 
Eq. (|J) is valid. The upper limit requires that the cutoff, 
r co ~ A/ _I/ , in the spatial extent of the avalanches is 
much less than the system size, or Af~ v <C L. Inserting 
Eq. (||) into this expression gives 

D = d+^—± . (5) 

The boundary case 7=1, realized for instance in the 
mean field version of the evolution model [^p| , ^llj5^ 1 , re- 
sults in an exponential relaxation to the critical attractor: 

A/ ~ e- As ' Ld , (6) 

where A is a numerical constant that is independent of L. 
This expression is valid as long as L d <C s <C L d InL, and 
the system reaches the stationary state when s ~ L d In L. 

It is straightforward to demonstrate, in a step by step 
fashion, that the gap equation (0) holds not only for 
the self-organization process in the evolution model, but 
also for the Sneppen interface model [5^| . For the self- 
organized LIM and the Zaitsev model, the distribution of 
internal forces in the system is not given by a flat distri- 
bution, unlike the evolution and Sneppen models. How- 
ever, the extremal dynamics allows one to define a fluc- 
tuating signal fmin(s) and therefore a gap function G(s) 
as above. Now, though, the average size in the jumps of 
the gap is not given simply by (1 — G(s))/L d . However, 
as long as the distribution of internal forces above the 
gap is not singular, when the avalanche is finished, then 
a weaker form of the gap equation holds, 

dG(s) 1 m 

ds L d <S > G{S) ' 1 1 

If the scaling ansatz Eq. (ra) holds, then depending on 
the value of 7, either Eq. (||) or Eq. (JsJ) describe the 
approach to the critical state. 

For d + 1 dimensional invasion percolation, the situa- 
tion is slightly more complicated. During the transient 
process, the size of the boundary where growth may occur 
is growing. It can be shown p5| that during the transient 
period of invasion percolation starting from a flat config- 
uration, the active boundary 6(s) of the invading cluster 
scales with the invaded volume s as 

b{s)/L d ~ (s/L d ) 9 = {s/L d ) i ^ . (8) 
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Here, ds is the fractal dimension of the boundary where 
growth may occur and D is the fractal dimension of the 
invaded cluster. The derivation of this equation is ex- 
plained in Appendix A. The system reaches a stationary 
state when s ~ L D . At this point, the size of the ac- 
tive boundary also reaches its finite size limit L dB . The 
proper gap equation for invasion percolation, taking into 
account the growth of the active boundary with time, can 
be written as 



dG(s) 
ds 



1 



L d {s/L d )a < S > G(S) 



(9) 



In the asymptotic critical region, the scaling ansatz Eq. 
(||) can be inserted into the gap equation (^). Integration 
gives 



A/ 



fc G( S ) = ( Td y 



(10) 



which holds for L d < s < L d+ ^J=ZJ . 

The scaling relations between these and other expo- 
nents are explained in the next section. One might 
naively suspect that the critical exponent 7 would be 
an independent exponent describing off-critical behavior. 
This is not the case. The results of the next section allow 
us to determine 7 in terms of the avalanche dimension D, 
and the avalanche distribution exponent, r, and solve the 
gap equation. 



IV. PROPERTIES OF THE STATIONARY, 
CRITICAL STATE 

In all previous SOC models, like for instance the BTW 
sandpile model jl6| or the earthquake models p8| , Self- 
Organized Criticality manifested itself in a power law 
distribution of bursts of activity, or avalanches, follow- 
ing a single perturbation. The new series of extremal 
SOC models is not an exception. In fact, avalanches 
in these models have an additional hierarchical struc- 
ture of subavalanches within avalanches. In this sec- 
tion, avalanches in the stationary state are defined and 
related to the avalanches in the transient process defined 
in the previous section. This enables us to determine 
the exponent 7, characterizing the transient, in terms 
of the stationary probability distribution of avalanches. 
The hierarchical avalanche structure is revealed in the 
"gamma" equation pqj . We then demonstrate the law 
for stationary states, 77 = ]39[ |. This law implies, 
among other things, that the fractal dimension of ac- 
tive sites is also determined by the stationary probabil- 
ity distribution for avalanches through the scaling rela- 
tion df = D(t — 1). All of the above mentioned results 
hold for forward avalanches; we close this section with 
a discussion of backward avalanches p6| , and a deriva- 
tion of a scaling relation for the Levy flight exponent, 
tt = 1 + Z>(2-t). 



Avalanches in the stationary state can be defined for 
any of the extremal models as follows: Suppose that at 
time s the smallest random number in the system was 
larger than f a , where < f a < 1 is an auxiliary param- 
eter used to define avalanches. Each of the new random 
numbers introduced at this time step can be smaller than 
f Q with probability f„ (for the flat distribution V{f)). 
According to the rules of the model, if one (or more) of 
the new random numbers is less than f Q , the smallest of 
those will be selected at the next time step (s + 1). This 
initiates a creation- annihilation branching process, where 
the sites with /, < f play the role of particles, and the 
particle with the smallest random number is chosen for 
updating at each time step. While the avalanche contin- 
ues to run, there is at least one such "particle" and the 
global, minimal random number is smaller than f a . The 
avalanche stops, say at time S + s, when the last "parti- 
cle" with a random number smaller than f a is eliminated. 
By definition, the global, minimal random number at this 
time step will be larger than f a once again. Thus, one can 
view the parameter f a as a branching probability for the 
creation of particles. We will call the avalanches, defined 
by this rule, / -avalanches. They can be easily extracted 
from the temporal signal of the model / m m(s), which is 
the value of the selected minimal number as a function of 
time s. In terms of this signal, the size of an / Q -avalanche 
is the number of time steps elapsed between subsequent 
punctuations of the barrier f a by the signal f m in ■ For the 
example given here, the size of the avalanche is S. The 
hierarchical structure of f a - avalanches is demonstrated 
in Fig. |. 

The statistics of the avalanches clearly depends on the 
value of the branching probability f . The larger it is, 
the larger is the expected size of the avalanche. In anal- 
ogy with ordinary percolation p3| , there must be a crit- 
ical value f c < 1 of the branching probability such that 
the expected size of the f c - avalanche cluster becomes 
infinite. That means that in the thermodynamic limit 
L — > 00 in the stationary state of the system, with prob- 
ability one, at least one "particle" with fi < f c will exist, 
and the signal f m i n will be smaller than f c , also with 
probability one. In analogy with ordinary percolation 
p3[ and other "branching processes" such as directed 
percolation |l4| , a variety of nonequilibrium lattice mod- 
els |64|] , or the propagation of an avalanche in the BTW 
sandpile model |ll|, we use the following scaling ansatz 
for the probability distribution P(S, /„) of /„- avalanches 
of size S: 



P(SJ ) = S- T g(S(f c - fo) 1 ^) 



(11) 



Here, r and a are model dependent exponents and g(x) 
is scaling function, which decays rapidly to zero when 
x 1 and goes to a constant when x — > 0. This scaling 
ansatz for various models has been confirmed by numer- 
ical simulations in @, @, @, (60), When 



8 



the auxiliary parameter f a is selected to be equal to f c , 
the avalanche distribution is a power law P(S) ~ S~ T 
without cutoff. When the parameter f„ is lowered below 
f c these critical avalanches are subdivided into smaller 
ones and acquire a cutoff S co = (f c — / ) -1 ' ff . The aver- 
age size of an / D -avalanche diverges as f a approaches f c 
as 

< S >~ (f c - /„)-* . (12) 

Simple integration of the power law (|ll]) gives the usual 
expression for 7 in terms of r and a, as occurs for example 
in ordinary percolation Q], fl63| : 

<S>= I SP{S, f a )dS ; 7 = — - . (13) 



A. The BS Branching Process. 



value of the gap, G, to the stationary f a — G- avalanche 
distribution. We emphasize, again, that the 7 for the 
stationary distribution in the evolution model, if it exists 
(i.e. if the scaling ansatz is valid), is the same 7 that en- 
ters into the solution of the gap equation for the transient 
behavior of the model. 

For the other models, we propose analogously that the 
G = f - avalanches during the transient also have the 
same scaling properties as the f a - avalanches in the sta- 
tionary state. The picture is that the scaling behavior 
of f a - avalanches is not affected by correlations at dis- 
tances larger than the correlation length £ ~ (f c — f )~ v 
set by f a . When the self-organizing system reaches a 
gap G = f , it has organized itself up to a scale £(/ ); 
at length scales smaller than this scale, it behaves as a 
critical system. Thus the exponent 7 that appears in 
Eqs. (|l2|, |l3]) is the same 7 that enters into the transient 
approach to the critical attractor, Eqs. (3-5). 



Unlike the other extremal models, the evolution model 
has an additional feature greatly simplifying numerical 
studies of f a - avalanches. The propagation of an Jo- 
avalanche in the evolution model is totally independent of 
the environment (the values of fa > /„), where it propa- 
gates. This means that, with respect to an f D - avalanche, 
all of the sites with /, > f a can be treated as vacuum 
sites. Sites are unimportant as long as they are not oc- 
cupied with the "particles" of the f a - avalanche. In order 
to study f a - avalanches we have to keep track only of sites 
that have random numbers fi < f - In simulations of this 
BS branching process p^J , which is exactly equivalent to 
an f a - avalanche in the evolution model, we first create 
2d+l random numbers, chosen from the flat distribution 
V, at the center of the lattice and its 2d nearest neigh- 
bors. Random numbers smaller than the parameter of 
the simulation, / Q , are stored along with their positions. 
At each time step, the minimal of the stored random 
numbers is "activated" until there are no more stored 
random numbers. At this point the f a - avalanche is fin- 
ished. This method gives an efficient way to study the f a - 
avalanche distribution and other properties, completely 
free from system size corrections. An avalanche is con- 
sidered infinite, and not included in the distribution, if 
its size is larger than a numerically imposed cutoff s max , 
which can be made arbitrarily large. Results for the ex- 
ponent t from simulations of the BS branching process 
are shown for one and two dimensions in Fig. ^. 

The equivalence of the off-critical BS branching pro- 
cess at a given value of f Q with an /„- avalanche in the 
stationary state of the evolution model also proves that 
f Q - avalanches in the stationary state of the evolution 
model have the same statistical properties as the G = f - 
avalanches during the transient; in particular they have 
the same probability distributions. The gap equation (|^) 
maps the transient, self-organizing behavior at a given 



B. The "Gamma" Equation 

We now proceed to establish a general relation for the 
number of sites covered by an avalanche as the critical 
state is approached. The relation is exact for the Sneppen 
and evolution models, and a similar exact relation holds 
for IP. It leads directly to a scaling relation between the 
exponents r and 7 valid for all extremal models. This 
scaling relation was previously derived for separate mod- 
els in g], ||, ||, H), & We wiU reproduce here in 
more detail the method of derivation used in p8| ], which 
gives not only the scaling relation but the exact values of 
coefficients in it. The following argument applies specifi- 
cally to the evolution and Sneppen models; we then gen- 
eralize it to the other models. 

Suppose that a value for the parameter f Q < f c has 
been chosen. The moments in time, s, when the global 
minimal number f m in{s) exceeds f a serve as breaking 
points dividing the temporal axis into a series of f - 
avalanches, following one after another. The probabil- 
ity that at any given moment the signal fmin(s) will be 
greater than f a is given by 

p(f min > fo) = l/(S) fo , for f <f c , (14) 

where (S) f a is the average size of an f a - avalanche. 

Let n cov denote the number of sites covered by an 
avalanche. These sites had their random number changed 
at least once during the course of the avalanche. Since 
each site can be updated more than once, n cov is a pri- 
ori different from the avalanche size S. In fact, for any 
avalanche 

n cov < AS , (15) 
where A is a constant which depends on dimension. 
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We can relate the divergences of these two physical 
quantities as the critical state is approached. This is ac- 
complished by noting that, at the moment in time when 
an f a - avalanche is finished, the random numbers at all 
n cov sites covered by the avalanche are uncorrelated and 
uniformly distributed in the interval [/ ,1]. This is be- 
cause at the last update for each site within the avalanche 
the random number was chosen from an uncorrelated, 
uniform distribution between [0,1] with the condition that 
it be larger than f a . We will make repeated use of this 
fact in what follows. 

Recall that when the f - avalanche started, the small- 
est random number in the system was larger than f . For 
simplicity, assume that this number was also larger than 
fo + df 0l where df will be taken to be vanishingly small. 
When the f a - avalanche finishes, each of the n cov sites 
has equal probability jzy- to fall within the interval df 

above f a . Therefore, n cov ^zj is the probability that at 
least one of the n cov sites left behind by the f Q - avalanche 
has a random number between f a and f a + df Q . If now 
the auxiliary parameter f a is raised by an infinitesimally 
small amount df , the breaking points that had f m i n (s) 
between f and f + df a will no longer stop the f + df - 
avalanches, and the average avalanche size will increase. 
This property of the avalanche hierarchy is demonstrated 
in Fig. 0. Consider a very long temporal sequence 
fmin(s) which is also a sequence of N, f a -avalanches. 
Increasing the auxiliary parameter to f a + df a , the num- 



ber N will be changed by dN 
to leading order in df for N — > 



-N < n c 



> 



dfo 



00. Of course, the sum 
of the temporal durations of these avalanches will remain 
constant. This leads to the following differential equation 
for the average size of an avalanche; 



d ^ n ( S )fa _ ( n COv)f 



dfo 



l-/o 



(16) 



This equation is exact for the evolution and Sneppen 
models in any dimension. It does not require the use 
of any scaling assumptions. 

In order to proceed further, we will now assume that 
the avalanche distribution obeys the scaling ansatz, Eq. 
(|i"T|). Then for f close to the critical value f c , the average 
size of the avalanche diverges as (/ c — / )~ 7 (Eq. ([jj])). 
Substituting this power law into Eq. (|l6|) gives 



7 =lim /„->/„ 



{n C ov)fAfc - fo) 

1-/0 



(17) 



If the critical exponent 7 exists, then Eq. ( |T7j ) is also 
exact. Note that the quantity 7 appears as a constant 
rather than a critical exponent. It represents the aver- 
age number of random numbers between f a and / c , left 
behind by an f a - avalanche that has died, in the limit 
that f — > f c . This number is universal. For example, 
it does not depend on whether one updates just nearest 
neighbors or also includes next nearest neighbors. 



A plot of (1 — f )/{n>cov) f , for different values of f Q , as 
shown in Fig. || for the two dimensional evolution model, 
gives f c as the intersection with the horizontal axis, and 
1 /j as the asymptotic slope close to f c . This enables us to 
determine f c very accurately, because the uncertainty in 
the measurement of (n cov ) and 7 leads to an uncertainty 
not in f c but in A/. Choosing a value of f a near the 
critical point, Af — f c — f is estimated from Eq. (p"^). 
Specifically, in the one dimensional evolution model we 
choose f = 0.665 and measure the average number of 
covered sites (n CO v)o.665 = 446.9 ±7.0. Using this method 
for the evolution model, we find f c = 0.66702 ± 0.00003 
for d = 1 and f c = 0.328855 ± 0.000004 for d = 2. Our 
one dimensional result agrees with Grassberger . As 
a by-product, the one dimensional result rules out the 
early speculation that f c is exactly 2/3 in one dimension. 
The exponent 7 can also be estimated accurately from 
points further away from f c . 

Figs. ^ and [l(] show the same plots for the one dimen- 
sional evolution model and the one dimensional Sneppen 
model. For the one dimensional evolution model, our re- 
sult, 7 = 2.70 ± 0.01, agrees with that of Jovanovic et al 
who measured 7 = 2.7 ± 0.1 E3| and Grassberger, who 
also measured 2.71 ± 0.03 @. 

For the one dimensional Sneppen model with L = 10 4 , 
we find fc = 0.465 and 7 ~ 2.13. This value differs 
substantially from that obtained by Sneppen 7 = 2.05 ± 
0.01 at f c = 0.46136±0.00005 at a system size L = 5xl0 5 
p7j. Our probably less accurate result is due to finite 
size effects. This value f c = 0.465 for the Sneppen model 
also does not agree with the value of f c = 0.4614 that we 
determine from the avalanche distribution (see Section 
IV (D)) for a larger system size L = 2 x 10 5 . 

For the evolution model, the value of f c obtained us- 
ing this new method, based on the "gamma" equation, 
is used as input to determine 7 from the divergence of 
the average avalanche size, as f a approaches f c , Eq. (p^). 
Our results for 7 for the evolution model using this con- 
ventional method (this was the method used in Refs. p2| 
and jy]]) are shown in Figs. |ll|, |l2|. The two different 
methods we used to measure 7 give consistent results for 
the evolution model. Having measured 7, we can calcu- 
late the exponent p for the relaxation to the critical state 
(Eq. (|)). We find p = 1/(7 - 1) = 0.588 ± 0.004 and 
p = 1.43 ±0.01 in the one and two dimensional evolution 
models, respectively. 

The scaling form, Eq. ([TT]) . for the divergence of the 
average avalanche size can be substituted into Eq. (pl|), 
which becomes p(f m in > fo) ~ (fc - fo) 1 , or alterna- 
tively, 



Pifn 



fo) ~ {fc fo)^ 



(18) 



From Eq. ( O ) it follows that the average number of 
sites {n CO v)f covered by an f Q avalanche scales near the 
critical point as 
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{n CO v)f ~ (fc - fo)~ 



(19) 



Since Eq. (^5|) holds for any avalanche, it is immedi- 
ately clear that 7 > 1 for any model which obeys Eq. 
(|l9|). This posteriorly confirms that Eqs. (||, ||) are the 
only allowed possibilities for the approach to the critical 
attractor via the gap equation. The case 7 < 1 is not 
physically realized by the models considered here. 

We are now in a position to consider the spatial ex- 
tent of avalanches. Since the interactions are via nearest 
neighbor, it is clear that all sites visited by an avalanche 
form a connected cluster. In one dimension any con- 
nected cluster is also compact. Therefore, in any one- 
dimensional model n cov ~ R cov , where R CO v is the length 
of the spatial interval formed by covered sites. 

In analogy with percolation, it is conventional to define 
the characteristic spatial size R of an avalanche cluster 
as the mean square root deviation of the set of all active 
sites in the avalanche from their center of mass. In this 
definition each site is counted with the weight given by 
the number of times it was visited by the avalanche. The 
avalanche mass dimension D is defined by the scaling 
relation 



R l 



(20) 



connecting the avalanche size S (temporal duration) 
to its spatial extent R. In the case of a compact d- 
dimensional set of covered sites, we can use another def- 
inition of the avalanche spatial size; 



R c 



(21) 



Assuming the absence of multifractal spatial scaling be- 
havior for avalanches implies that 



Rc 



R 



(22) 



In all extremal models that we have studied numer- 
ically thus far (excluding IP which is discussed at the 
end of this section), the avalanche mass dimension D 
was measured to be larger than the dimension of space 
d, and the spatial projection of all active sites within the 
avalanche was observed to form a compact object of di- 
mensionality d. In what follows, we assume this is true 
for every extremal model below its upper critical dimen- 
sion. Above the upper critical dimension (provided that 
it exists) the fractal dimension of the collection of cov- 
ered sites is given by its mean field value d uc , and can 
no longer form a compact object in the space of dimen- 
sionality d > d uc . Thus in all of the scaling relations that 
include dimensionality, derived below and summarized in 
Table GL d should be replaced with d uc for d > d uc . This 
behavior is somewhat analogous to hyperscaling relations 
in equilibrium critical phenomena. 

In the case compact avalanches | p8[ , n cov = R d ov ~ 
Rd ^ s d ' D and from Eq. @ it follows that (A/)" 1 ~ 



result, 



) ~ jf' ^ S d ' D S- T dS ~ ^{r-dlD-l)l^ Ag a 
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a = 1 + d/D - t 
(2-r) 
(1 + d/D - t) 



(23) 
(24) 



The spatial correlation length exponent v describes 
the scaling of the cutoff in the spatial extent of the 
avalanches, r co ~ A/ - ". Since s co ~ r® , 

v = l/aD= l/(D + d- Dt) . (25) 



It is interesting to observe that Eq. (|23j-|25|) imply that 
the self-organization time to reach the critical attractor 
is independent of the initial state of the system, for al- 
most any initial condition. A system of size L reaches 
the stationary state when (A/(s))~" ~ L. The time s org 
required for this scales as s org ~ L D , where D = d+ 
(Eq. H). Substitution of Eqs. ( p3[j25| ) into this expression 
for D gives the very simple result s org = L D . or D = D. 
Instead of starting with the initial condition where the 
random numbers are uniformly distributed in the inter- 
val [0,1], as assumed in the gap equation (Eq. ||), we may 
start the organization process in a state where all /'s are 
set equal to 1 except one site which has the value of / 
equal to, say, 0.99. In this case, none of the original l's 
are ever selected as the minimal site. The organization 
process is finished at the moment when the last 1 is de- 
stroyed. This clearly takes ~ L D time steps, the same as 
for the transient that is governed by the gap equation. 

We checked these scaling relations numerically. Figs. 
|l3|a, |l4| a show measurements of the avalanche dimension 
D for the one and two dimensional evolution model. At 
the critical point f c , R cov was measured as a function 
of S for running avalanches and their subavalanches in 
the BS branching process, and averaged over more than 
10 11 updates. In d = 1, D = 2.43 ± 0.01, and in d = 2, 
D = 2.92 ± 0.02. These values are consistent with the 
scaling relation Eq. (^i|) and our measured values for 
t and 7. In fact, Jovanovic et al [fi"2"f measured the ex- 
ponent for the divergence of the average spatial size of 
the avalanche to be 0.98 ± 0.03 (i>± in their notation) 
for the one dimensional evolution model. Taking into 
account that in one dimension the spatial size is simply 
proportional to n cov this result is consistent with the ex- 
act value of 1 we predict. In one dimension, Eq. ( p5|) can 
be written as 

r - 1 = av - a . (26) 

For the one dimensional evolution model, the results from 
Ref. @],r= 1.08±0.05,cr = 0.35±0.02andl/D = cri/ = 
0.43±0.01 are consistent with this last relation. Similarly, 
Grassberger's results r = 1.073±0.003, a = 0.343±0.004 
and l/D = av = 0.4114 ± 0.002 Q are also consistent 
with Eq. (p6|) - 
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This agreement occurs in spite of the fact that it is 
very difficult to determine the true asymptotic value of 
D = X/(av). One reason is that the apparent exponent 
may be nonmonotonic, as demonstrated in Fig. p"3|b. The 
effective exponent as plotted here is measured by numeri- 
cally calculating the local slope of the curve R vs. S. It is 
apparent for the two dimensional model that the effective 
exponent, av = 1/D, does not approach its asymptotic 
value monotonically but first overshoots with a maximum 
around S ~ 10 and then undershoots with a minimum 
around S ~ 10 2 . There are indications that similar be- 
havior may occur in the one dimensional model at much 
larger avalanche scales. Our numerical capability for the 
one dimensional model is not sufficient to be sure that 
we are observing asymptotic behavior despite studying 
3 x 10 11 time steps, and thus 3 x 10 11 avalanches and 
sub- avalanches. In addition, simulations of the related 
M —* oo evolution model |53[| , where D = 4 is known ex- 
actly, indicate that there are logarithmic corrections for 
Rcov Thus, the numerical studies thus far have not con- 
clusively refuted the conjecture fs^ , Q, based on sym- 
metry arguments, that the exponents D and r in the evo- 
lution are the same as directed percolation, or Reggeon 
Field Theory. It seems dangerous to us to base conclu- 
sions about universality on differences of less than two 
percent in the observed "asymptotic" regime at S ~ 10 6 
pl[ in one dimension. In two dimensions, the disagree- 
ment is within error bars and less than one percent. 

Eq. ( |l6| ) can be generalized to include IP with the fol- 
lowing consideration. In IP the list of growth sites from 
which the minimal number is selected is not a compact 
e?-dimensional lattice of linear size L, but the fluctuat- 
ing fractal boundary of the invaded cluster. Comparing 
the list of potential growth sites before and after an Jo- 
avalanche, one observes that some sites were invaded dur- 
ing this specific avalanche and thereby removed from the 
list. In addition, the result of an / -avalanche in IP-2 
and IP-3 may be a new internal "lake" in the invading 
cluster, so sites on the lake perimeter are now excluded 
from the active boundary as well (see Appendix). On 
the other hand, some new boundary sites may have been 
added. We will denote this latter number of new bound- 
ary sites added by an avalanche of size R as n new (R). 
Since these sites constitute a part of the active boundary 
of the region invaded by the avalanche, it is reasonable 
to assume that 



,(R) ~ R d 



(27) 



where (is is the fractal dimension of active boundary of 
the invaded region. All of these new random numbers are 
uncorrelated and uniformly distributed between f a and 
1. By replacing the number of covered sites, n cov , with 
the number of new sites added to the boundary, n new , 
we find that Eqs. (]l^-pT|) apply to IP as well. 

Using Eq. (p7|) we can derive exponent relations anal- 



J*f "° s d B /D s -r ds _ A j(r-l-d B /D)uD^ wh[ch giveg 

d B - l/v 



T = 1 + 



D 



(28) 



(29) 



Although the exact equations (|16|-|l7D, with n cov replaced 
with n new , are new for invasion percolation, the scal- 



7= 1 + (D- d B )v 



ing relations (|2J|-|29|) are well known. To the best of our 
knowledge, they were first derived by Gouyet [f39| . Sub- 
stitution of Eqs. ( p9|) and (||) into Eq. (10) simplifies the 
final form of the transient approach to the critical state 
for IP: 



( S /L d y 



•HD-d B ) 



Af ~ (s/L a ) 



(s/L d Y 



d\-i 



(30) 



The time required to complete self-organization s c 



L d+ ^T) ~ L D . Eq. (|0|) agrees with results from ||] 
obtained using different arguments. 

Eq. ( ^4]) holds for any extremal model where f m in(s) 
has a well defined threshold level / c in the thermody- 
namic limit. At the same time, though, Eqs. ( [l(^|l7j ) are 
not valid for the Zaitsev and LIM models. The reason 
is that for these two models the distribution of internal 
forces left by an f a - avalanche is not flat and uniform in 
the interval [/ , 1]. However, as long as this distribution 
is not singular in the limit that /„ — > / c , and these inter- 
nal forces are uncorrelated, as is reasonable to assume, 
the result (n cov ) ~ (/ c — / ) _1 still holds, and, therefore, 
the scaling relations (|23]-p5|) are valid. For interface mod- 
els, it is customary to express scaling relations in terms 
of D, v as basic exponents. Our scaling relations will 
then read 



T = 1 



1/V 



D 



7 = 1 + (D - d)v 



(31) 
(32) 



ogous to Eqs. <Mm for IP: (A/)" 



>>/„ 



C. Law for Stationary States: 77 = 

In this section we derive a law for the activity in the 
stationary state for the class of SOC extremal models 
defined in Section II. Recall that an f c - avalanche starts 
when all sites in the system have random numbers, or 
internal forces, above the threshold f c . In the follow- 
ing discussion, we will only consider f c - avalanches, and 
therefore drop the label f c . During the course of the 
avalanche there will be sites, denoted as active sites, or 
"particles" , where the random numbers fi are less than 
the threshold f c . The number of these active sites after 
s updates is denoted the "activity", n(s). This quantity 
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is just the number of sites in the system that have ran- 
dom numbers less than f c . The current avalanche stops 
when all sites in the system again have random numbers 
above f c , or when n(s) first returns to zero. Then a new 
avalanche will start somewhere else in the system. The 
probability distribution for avalanche sizes, P(S) ~ S~ T , 
is therefore, by definition, also the probability distribu- 
tion of interval lengths between subsequent returns of 
n(s) to zero. 

Fig. |l6| shows a typical n(s). The collection of return 
points, {n(s) = 0}, forms a fractal with dimension d — 
t — 1 on the one-dimensional time line. This relationship 
for d is explained in the next section (see Eq. (|5(J)). As 
long as d < 1, the threshold f c is well defined in the 
thermodynamic L — > oo limit. The initiation of a new 
avalanche in the stationary state can be viewed as the 
injection of a single "particle" into the system. Thus, the 
average number of return points, or number of injected 
particles, added in an interval of s steps is nj^j ~ s d . 

We examine the average number of active sites, S steps 
after a single particle has been injected into the system, 

< n(S) >~ S v . This quantity is by definition the av- 
erage activity tisurv(S) ~ S ds of the avalanches that 
survive S steps multiplied by the probability of survival, 
P(s' > S) ~ S 1 ~ T . Avalanches that die before S steps 
are counted with zero particles in the average < n(S) >, 
while the surviving avalanches are counted with their ac- 
tual number of active sites n(S). The quantity nsuRv{S) 
is an average over avalanches that survive S steps, while 

< n(S) > is an average over all avalanches, hence the 
factor P(s' > S). We now utilize the hierarchical nature 
of the stationary process: we propose that the activity of 
the large avalanches after S steps scales in the same way 
as the amount of activity injected into the system during 
a time interval of length S: 

nsuRv(S) ~ n INJ (S) , or d s = d = t - 1 . (33) 

That is, there is only one dimension for the activity in 
the critical process. Thus, 

< n(S) > = nsuRv(S) P(s' > S) 
~ n INJ {S)S 1 ~ T ~ S° ; i.e. 77 = . (34) 

In order to illustrate the argument by an elementary 
example, consider a simple, uncorrelated one dimensional 
random walk, with the activity increasing or decreas- 
ing by one with equal probability at each step. This 
is the mean field, SOC state in many models pCj], [ pjl , 
|p2| . The infinite random walk can be viewed as a se- 
quence of avalanches. These avalanches are, by defi- 
nition, the intervals between subsequent returns of the 
walk to zero. The exponent r for the first return times 
is t = 3/2, so the probability of s' exceeding S scales 
as S^ 1 ' 2 . The average number of returns in S steps, 
niNj(S) ~ S 1 / 2 , scales in the same way as the activity 



of surviving avalanches after S steps, nsuRv(S) ~ S 1 / 2 . 
Thus, the average activity after S steps of a single 
avalanche is < n(S) >~ S 1 I 2 S- 1 ' 2 ~ 5°, and rj = 
rigorously. 

Fig. [I?] shows the average activity of a single 
avalanche, < n(S) >, in the two dimensional evolution 
model. The exponent r\ reflects the average activity S 
time steps after starting a single realization of the BS 
branching process at f c . The quantity < n(S) > sat- 
urates for large S, varying only 2% over four decades: 
this is consistent with rj = ± 0.002 in two dimensions. 
Fig. |l^ is our strongest numerical confirmation of the 
rj = prediction. As eluded to earlier, c.f. Fig. [l4|b, in 
one dimension there are larger corrections to scaling and 
the saturation regime obtained numerically is smaller, al- 
though the overall behavior is the same, as shown in Fig. 
|l8| . The slope is decreasing and reaches n = 0.01 ± 0.02 
at s = 6 x 10 7 . 

The fractal dimension of active sites within a surviving 
avalanche is defined by the relation nsuRv( r ) ~ rd - f > 
where r is the spatial extension of the avalanche cluster. 
Using the relation s ~ r D , and Eq. (|33]), we find the 
relation 

d f = Dd s = D(t - 1) , (35) 

connecting the geometrical properties of the active sites 
to the properties of the entire avalanche clusters. Fig. |l7| 
(|l8|) also shows nsuRv(S) vs. S in two (one) dimensions 
for the evolution model. At the largest time scales mea- 
sured, the slope of the curves yields d s = 0.11 ± 0.02 and 
d s = 0.25 ± 0.005 in one and two dimensions. In two 
dimension, this is consistent with our measured value 
of r - 1 (c.f. Fig. but in one dimension our mea- 
sured value of t — 1 is somewhat lower. This is proba- 
bly a consequence of the slow convergence of the slope, 
as is evident from Fig. |l^. Actually, Grassberger mea- 
sured nsuRv(S) ~ S° for the one dimensional evolu- 
tion model, which agrees with his measured value for r 
plf, and thus with Eq. @ and <q = 0. 

The usual dynamical exponent relating space and time, 
z, where t ~ r z ~ s z / D , can also be measured. The par- 
allel time unit, t, is different than the sequential time 
unit s. It is defined as the average number of update 
steps for changing all active sites, so that time increases 
by an amount l/n(s) at each update; i.e. s — > s + 1, 
t — ► t+ 1/nif). It appears naturally in models for depin- 
ning at constant force, for example, because all unstable 
sites move together. We introduce the dynamical expo- 
nent z also for the sequential models in order to compare 
dynamical critical behavior in the different situations. 

For this definition of time, clearly s ~ nt and d s + 
z/D = 1. Within the superuniversal class where 77 = 0, 
including the extremal SOC models discussed here, the 
exponent z is given by the scaling relation 

z = D — df = D(2 — t) . (36) 
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We measured the dynamical exponent z — 1.19 ± 0.05 
for the Sneppen model in one dimension as shown in Fig. 
|l9| . This value is close to the theoretically predicted value 
z = 1.21, obtained by using Eq. @ with D = 1.63 ||, 
and our measured value of r. We have also measured 
the dynamical exponent z for the one and two dimen- 
sional evolution models as shown in Fig. |2^. In two 
dimensions, z = 2.16 ± 0.02. Inserting our measured 
values for r and D into the expression for z we predict 
z = 2.20 ± 0.04, in fairly good agreement. In one dimen- 
sion we hnd z = 2.26 ± 0.05 to be compared with the 
somewhat lower measured value 2.10 ± 0.05. This dis- 
crepancy again may reflect slow convergence for the one 
dimensional evolution model. Jovanovic et al ft2| | mea- 
sured a dynamical exponent (in their notation v\\jv±) 
= 2.1 ± 0.1 for the one dimensional evolution model us- 
ing the tree structure of the avalanches. The dynamical 
exponent that they measured may be equivalent to our 
definition of z, and their result is consistent with our 
prediction, again within numerical uncertainty. 

An alternative way of deriving r; = follows: Consider 
either the evolution or Sneppen model in the station- 
ary state, and select an arbitrary value J below J c . We 
consider the number of active sites with random num- 
bers below J c . As explained in the preceding subsection, 
when an J a - avalanche ends, it leaves behind uncorre- 
cted random numbers uniformly distributed between /„ 
and one. As a result, the average number of active sites 
created, or left behind, by an f a - avalanche is exactly 
< n cov > fa (f e - / )/(l - J ) (see Eq. (17)). We as- 
sume that when the J a - avalanche started, the active sites 
between J a and J c were distributed on a fractal, with 
dimension df. As a result, the average number of ac- 
tives sites that are destroyed by the J a - avalanche scales 
as < R di > f o , where R is the spatial extension of the 
avalanche. In the stationary, critical state the average 
number of active sites created and destroyed by an J a - 
avalanche must be equal. Assuming a power law distri- 
bution of avalanches sizes with cutoff, R co ~ (J c — J ) v ', 
we find 

l=v{d-d f ). (37) 

Using the scaling relation, r = 1 + (d — l/v)/D (Eq. 
(p5|)) we again obtain our fundamental result r\ — df /D — 
t + 1 = 0. Eq. (|37|) is reminiscent of the hyperscaling 
relation [3 = v(d— df) valid for equilibrium systems, and 
for directed percolation |]70| . 

In order to take this alternative argument over to the 
Zaitsev model and LIM, one makes the same assumption 
that was made for the gamma equation, i.e. the distribu- 
tion of internal forces left behind by an J Q - avalanche is 
uncorrelated and smooth. We have measured < n(s) > 
and nsuRv(s) for the Zaitsev flux creep model p9[ in 
one dimension, as shown in Fig. These results are 
consistent with t] — 0. Finally, we predict that 77 = for 
invasion percolation. 



Although the 77 = law applies to a broad class of 
SOC phenomena, it does not generally apply to systems 
that are tuned to be critical, since for non-SOC systems 
the critical state may not be stationary. For instance, 
in directed percolation 77 = df / D — r + 1 ~ .21 |nj in 
one dimension. Also, in directed percolation the fractal 
dimension of active sites obeys a hyperscaling relation, 
df = d — D(t — 1). This hyperscaling relation is incon- 
sistent with the relation df = D(t — 1) that holds for 
the extremal SOC models that we study. Thus directed 
percolation does not belong to the superuniversal class 
defined by the dimension independent law rj = 0. 

The models for depinning at constant force, such as 
TLB or parallel LIM, also do not in general obey r] = 0, 
or the corollary results df = D(t — 1), z = D(2 — r). For 
example, in the one dimensional TLB at constant force, it 
has been predicted and confirmed numerically that z — 1 
p2| . Our measured value for this model z = 1.00 ± 0.05, 
as shown in Fig. [l^, is consistent with this prediction. 
Since D ~ 1.63 and r ~ 1.26 are the same as for the 
Sneppen model (see section VI), one can compare with 
the different prediction z ~ 1.21 based on the 7; = 
law. Since 77 is not zero for the tuned depinning models 
at constant force, while it is zero for the SOC depinning 
models, the critical dynamics of a self-organizing system 
are different than the critical dynamics of a tuned sys- 
tem. Thus SOC cannot be simply viewed as sweeping an 
instability . If the TLB model were to be studied at 
constant velocity, we predict that 77 = and z ~ 1.21, as 
for the SOC (Sneppen) case. Thus the dynamical criti- 
cal exponents, such as z, depend on how the depinning 
transition is probed. 

The existence of a stationary limit may imply that 
r/ = for a very broad class of SOC models, beyond 
those studied here. In addition, Pietronero and co- 
workers explicitly make use of the stationary limit in the 
Fixed Scale Transformation method. This suggests that 
stationarity may potentially provide more powerful tools 
to understand SOC and fractal growth phenomena in a 
wider range of systems. 

D. Backward Avalanches. 

The dynamics of extremal models forms a hierarchical 
structure. We have defined J a - avalanches as the activ- 
ity between subsequent moments in time when the signal 
Jmin{s) is larger than an auxiliary parameter J Q . Then, 
as the parameter J a is decreased, larger avalanches are 
subdivided into smaller ones because new breaking points 
appear. For the sake of clarity in this section, we will 
now refer to avalanches defined by this specific rule as 
/ -punctuating avalanches, since they correspond to a 
sequence of events between subsequent punctuations of 
the / -barrier by the signal J m i n (s) • Another way of look- 
ing at the hierarchy is to define avalanches by a slightly 
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different rule. According to the new rule, a forward 
avalanche begins at every time step s. It runs as long as 
f minis + s') < f minis), and will stop at the first moment 
s + S forward in time when f minis + S) > f minis)- This 
avalanche is similar to an Jo-punctuating avalanche, ex- 
cept that now we set f a exactly to the value of the signal 
f minis) at the starting point of the avalanche. In the evo- 
lution model, the mapping to the BS branching process 
proves that the probability to have an / Q -punctuating 
avalanche of size S, is exactly the same as the condi- 
tional probability to have a forward avalanche of this size, 
given that the value of the signal at the beginning of the 
avalanche was f . That is, the exact knowledge of the 
value of the signal f minis) > fo at the starting point of 
an / -avalanche doesn't influence its statistics. This fol- 
lows from the observation that the site carrying f minis) 
gets a new random number at the very first time step of 
the avalanche and all the information about its previous 
value is erased from the system. Therefore, the condi- 
tional probability Pf(S,f ) to have a forward avalanche 
of size S, given that the signal at its starting point was 
equal to f , is exactly equal to the probability P(S, f ) to 
have an / G -punctuating avalanche of size S. Then from 
Eq. «, 



PfiS, f D ) = P(S, f ) = S- T giSif c - f 



(38) 



For the other models we lack a rigorous proof, so we 
simply conjecture that both PfiS, f a ) and P(S, f ) scale 
(e.g. Eq. y) with the same exponents r and a, but with 
possibly different scaling functions gfix) and gix). With 
this assumption, the results which follow in this section 
also apply to the other SOC extremal models. 

To study the properties of the signal under time 
reversal it is useful to define backward avalanches. 
Now we look for the first moment back in time when 
f mmis - S) > f minis) = fo- The definitions of both for- 
ward and backward avalanches are illustrated in Fig. |2^. 
This figure demonstrates the hierarchy in the avalanche 
structure: all forward and backward avalanches that 
start inside one big forward (backward) avalanche are 
constrained to not go beyond the limits of the parental 
avalanche and, therefore, can be considered to be its sub- 
avalanches. Each sub-avalanche in turn has its own sub- 
avalanches, and so on. One can look at the entire activ- 
ity in an extremal model as one great parental critical 
avalanche, which began in the distant past. It contains 
sub-avalanches of all sizes. 

Since all extremal processes are intrinsically irre- 
versible, it is possible to have different statistical proper- 
ties for the forward and backward avalanches. In analogy 
with Eq. (H|), we conjecture a scaling form for the con- 
ditional backward avalanche probability distribution; 



Pb(SJo) = ^S- T »g b (S(f e - fo) 1 '") 



(39) 



zero for x ^> 1. We will prove later that the cut-off ex- 
ponent a in this expression is the same for both forward 
and backward avalanches, while the power law exponents 
Tf, and r are different. In fact t < 1, so that a normal- 
ization factor l/N must be included in the conditional 
probability distribution for backward avalanches. 

According to Eq. (W), the minimal numbers selected 
at different time steps are distributed with probability 
density p(f min = f D ) ~ ifc - fo) 1 ' 1 , where 7 = { s 
the critical exponent governing the divergence of the av- 
erage punctuating avalanche size. The distribution of 
all forward avalanches, P^ ll iS) is obtained by integrat- 
ing the conditional probability from Eq. ( p8[ ) with the 
proper weight from Eq. ( |l8| ) to give 

ffc 

Pf iS) = / Pf iS, f ) Pifmin = fo)dfo 
JO 

/C S- T giSif c - f ) 1/a ) if c - fo)^ 1 df 



The unexpected result is that the exponent for this distri- 
bution has the superuniversal value -2 in all dimensions. 

For the distribution of all backward avalanches t < 1 , 
and the normalization factor l/N — 1/ S~ rb 5b(5 , (/ c — 

fo) 1/,T )dS = if c - fo)^ enters. Integrating Eq. (|3~ 
gives 



P b all iS) = J PbiSJo) pif mm = f ) df 
^ g—r b —(r(l—T b )/a—cr(2—T)/a _ qr—3 



(41) 



where gbix) is a scaling function that rapidly decays to 



The exponent t does not enter into the final expression 
for PfriS) as long as Tf, < 1. The exponent for the dis- 
tribution of all backward avalanches is model and dimen- 
sion dependent and is related to the usual f a - avalanche 
distribution exponent t. 

In numerical simulations, the exponent r is conven- 
tionally measured from the distribution P(S*, f ) of Jo- 
punctuating avalanches with the value of f carefully cho- 
sen as close as possible to the actual threshold f c - Since 
backward avalanches are defined at every time step, while 
PiS,f ) gets contributions only when f m i n (s) > f a , 
the distribution Pj? ll (S) has much better statistics than 
PiS, fo) after the same number of time steps. It is also 
very convenient that P^ ll iS) automatically has no cut- 
off, so one does not need to know f c in order to measure 
its power law. Thus, P^iS) provides, in principle, a 
very accurate way to measure r for the extremal SOC 
models discussed here. However, we have not yet pushed 
this technique to its ultimate limit. 

We have measured t£ u = 3 — r in the one dimen- 
sional Sneppen model, and the one and two dimensional 
LIM. These results are shown in Figs. ^3] and ^4|. For the 
Sneppen model, our numerical results from the backward 
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avalanche distribution give r = 1.255 ± 0.02, in agree- 
ment with the value of 1.25 ± 0.05 measured by Tang 
and Lcshhorn ]54j] , and r = 1.26 ± 0.01 measured by us 
(Fig. ^5|) using the conventional method. The LIM re- 
sults, T XD = 1.13 ± 0.03 and t 2D = 1-29 ± 0.03, are in 
agreement with the results from numerical simulations 
on self-afhne one dimensional interfaces in porous media 
from fl6Cfl , and two dimensional domain walls in the Ran- 
dom Field Ising Model Q . This supports the possibility 
that these models might be in the universality class of the 
LIM; a more detailed comparison is made in Section VI. 

A study of backward avalanches leads to an exponent 
relation for r b , and a proof that / Q -backward avalanches 
have the same cutoff as / Q -punctuating (or forward) 
avalanches. For the sake of simplicity we concentrate 
on the case of the evolution model, where the only as- 
sumption is that the scaling forms, Eqs. (38, [39]) for the 
avalanche distributions exist. With the additional as- 
sumptions that have already been mentioned, the result- 
ing scaling relations also apply to the other SOC extremal 
models defined in Section II. 

Consider an arbitrary / Q -punctuating avalanche. The 
probability distribution P(S, f ) for the size S of this 
avalanche is given by Eq. (^) . For this f a - punctuating 
avalanche to be also a valid / -backward avalanche, start- 
ing at s + S and running backwards in time down to s, 
one needs fminis + S) = f a . We next calculate what frac- 
tion of / -punctuating avalanches are also / -backward 
avalanches. Suppose we have a temporal sequence 
f min{s) which is an ensemble of N, / D -punctuating 
avalanches, where iV is a sufficiently large number. The 
average number of f a - punctuating avalanches of size S 
in such an ensemble is given by N(S) = NP(S,f ). At 
the end of any / -avalanche of size S, n cov sites have ac- 
quired new random numbers. If the avalanches are com- 
pact, then n cov ~ S d ^ D . All these random numbers are 
uncorrelated and uniformly distributed between f a and 1. 
To have f m in(s + 5) = f we need the minimal number 
in the system to lie between f Q and f + df . This num- 
ber can be only at one of these n cov updated sites, since 
at the beginning of the avalanche every number in the 
system was larger than f a . The probability that at least 
one of these numbers will be between f Q and f a + df is 



given by n c 



dfa 



gd/D 



dfa 



The number N b {S) of 



v cov X—f ' ~ ^ 1 — / 

valid /o-backward avalanches of size S in our ensemble is 
N b (S)dfo = NP(S,f )n cov ^- o ~ NP(SJ )S d / D ^f- o . 
Therefore, the conditional probability distribution of f a - 
backward avalanches obeys: 



P b (S,f )~S d / D P(SJ ) 



(42) 



where the proportionality constant is determined from 
the normalization condition X)s=i Pb(S, f a ) = 1. 

Intuitively, it is clear that larger / Q -punctuating 
avalanches affect a larger number of sites. These larger 
avalanches are more likely to leave behind a new ran- 



dom number between f a and /„ + df and thus consti- 
tute a valid / -backward avalanche. This explains why 
the distribution of backward avalanches acquires an ad- 
ditional factor of S d / D compared to punctuating (or for- 
ward) avalanches. Eq. ( f42| ) shows that both forward and 
backward avalanche distributions have the same cut-off 
as a function of / D , and their power law exponents obey 
the relation r b = r — dj D. 

FromEq. (gl|), we get a particularly simple expression 
for r b : 

r b = l + -^-- = l-- = l-. . (43) 



Since a must be positive, this proves that r b < 1 and 
posteriorly confirms the validity of Eq. ( j4l"| ) . 

All equations in this section apply to IP with the di- 
mension d replaced with ds (see Appendix A). The im- 
portance of forward avalanches in IP was first recognized 
by Roux and Gouyon |35[ |. However, they made some 
erroneous assumptions which lead to incorrect scaling re- 
lations. Our exponent relations agree well with their nu- 
merical results, though. For IP-2, where ds = 1-75 and 
D = 1.89, they measured r b a " = 1.50 ± 0.04 which is 
consistent with our prediction t£ u = 1.47 based on Eq. 
( [43| ) and the assumption that v = 4/3 as in ordinary 
percolation [|63[. 



E. Levy Flight Distribution 

In the stationary state, the minimal site jumps 
throughout the system in a correlated and anomalous 
fashion which has some similarity to the usual Levy 
flight picture. Specifically, one can record the spatial 
location f m i n (s) of the current extremal site (with the 
smallest random number) as a function of time s [p9[ , 
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The distribution Pj um p(r) of jumps 



\Tmin (s) — f m in {s ~ 1) | between subsequent extremal sites 
follows a power law: 



P 



jump 



('•) 



(44) 



This behavior is reminiscent of the Levy Flight Random 
Walk (LFRW), where at every time step the walker jumps 
in a random direction by a distance that is drawn from 
a power law distribution. In contrast to the uncorre- 
lated LFRW, the jumps of activity in extremal models 
are temporally correlated, so that even if n > 3 and, 
therefore, the jump distribution has finite first and sec- 
ond moments, the process does not necessarily reduce to 
ordinary diffusion. 

The exponent tt can be related to other exponents as 
follows: Consider a backward avalanche that started at 
time step s. Suppose it has size S. By definition, f m i n {s— 
S) > f minis), while f m in(s -k) < f minis) for 1 < fc < 
5 — 1. Looking at the same sequence of events forward 
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in time from time step s — S to s, one notices that the 
forward avalanche with f a = f m i n (s — S) that was started 
at time step s — S is still running at s. At this moment 
s, n cov ~ S d t D sites acquired new random numbers since 
time step s — S. The active sites at time steps s and s — 1 
are both rigorously constrained to belong to this specific 
set of covered sites. 

Looking further back in time one may find another 
forward avalanche, which is still running at time step 
s. Say, it started at time step s — S' < s — S with 
fmin(s ~ S') > fmin(s — S), and has covered a larger spa- 
tial region R' ~ S n ^ D . Such an avalanche will contain 
the avalanche that started at time step s — S as one of its 
sub- avalanches. Active sites at time steps s and s — 1 will 
obviously belong to the bigger spatial region of size R' as 
well. The importance of the backward avalanche is that it 
automatically selects the smallest forward avalanche con- 
taining both sites r m in(s — 1) and r m i n (s) and, therefore, 
imposes the most restrictive constraint on their relative 
positions. 

For the evolution model, we proceed by showing that: 
1) the position of activity f m i n {s) at time step s is equally 
likely to be at one of the n cov sites with new random num- 
bers; 2) It is uncorrelated with the previous position of 
the active site r m j„(s-l), within the set of n cov sites. To 
do this we will again use the powerful observation that 
any / Q -avalanche of size S leaves in its wake a set of un- 
correlated random numbers. For our avalanche in ques- 
tion we can set f Q equal to f m in{s)- Then n cov new num- 
bers are uniformly distributed between /„ and 1 and are 
equally likely to host the current global minimum. Look- 
ing at the n cov sites at time step 8 — 1, just before the up- 
date was performed, one observes that almost all of these 
sites have already acquired their final uncorrelated num- 
ber between f a and 1. The only sites which can poten- 
tially be active and, therefore, correlated are r m j n (s — 1) 
itself and its nearest neighbors. Because at the next time 
step the f Q - avalanche dies out, all 2d + 1 new random 
numbers created at this last time step must be larger than 
f a . They simply join the rest of the n cov sites, which al- 
ready have their random numbers between f a and 1, and 
become indistinguishable from them. Therefore, the po- 
sition of the active site at time step s is uncorrelated from 
the particular position, within the set of n cov sites, of the 
active site at time step s — 1 . This finishes the proof that 
the current jump of activity r — \r m m(s) ~ r m in(s — 1)| 
for the evolution model is bounded only by the linear 
size R of the backward avalanche that starts at time step 
s. We propose that this bound also holds for the other 
extremal models considered here. 

Given this bound, the probability Pj U mp( r > Ro) to 
have a jump distance r larger than R a for large R scales 
in the same way as the probability to have a backward 
avalanche of linear size R larger than R a , or, alterna- 
tively, volume S larger than R®. Therefore, Pjumpif > 



R ) ~ Pg u {S > R D ) = R D(lb L >. Substituting the 
expression for t£ u from Eq. ( |4l| ) into this relation and 
differentiating with respect to R a , we get the final ex- 
pression for the exponent tt, 



TT = 1 + D(2 - t) 



(45) 



This expression is in agreement with the result tt = 1 + 
r f/v= 1 + vD{2 — t)/v = 1 + D{2 — t) which was derived 
for the Sneppen model using different methods in p5[ . 
For invasion percolation, the relation between tt and rf 11 
was first derived by Roux and Guyon |39] . Unfortunately 
they had an erroneous expression for t^ 11 . 

Based on Eq. (0) and our numerical values for D and 
r, we predict 7r = 3.26 (n = 3.20) in the one (two) dimen- 
sional evolution model. We have measured the exponent 
7T = 3.23 ± 0.02 (Fig. |(|, and Jovanovic et al ff| also 
measured 7r = 3.1±0.2in the one dimensional evolution 
model. In the one (two) dimensional LIM, it = 2.93±0.01 
(Fig. |27|) ( TT = 2.89 ± 0.03 (Fig. |§)). These results 
also agree with Eq. ( f45| ) and the numerical values we 
obtained for the exponents D and r. For the one dimen- 
sional Sneppen model, we predict it = 2.21; Sneppen and 
Jensen measured tt = 2.25 ± 0.05 ]72] |, while Tang and 
Leschhorn measured tt = 2.20 ± 0.05. In two dimensions, 
Falk et al [j73j measured n = 2.2 ±0.2 consistent with our 
prediction based on their measured values x — 0.50±0.03 
and t = 1.45 ± 0.03. For IP-3 Furuberg et al measured 
tt = 2.1 ±0.1 [Q, which is consistent with their measured 
value d B = 1.37 and D = 1.82 and Eqs. @ and @. 



V. THE FRACTAL PATTERN OF ACTIVITY 

For the SOC extremal models discussed in this article, 
the dynamics consists of a series of extremal events fol- 
lowing one after another. At any given time step, s, there 
is one and only one lattice site where activity occurs. The 
extremal character of this dynamics lies in the fact that 
this site is always selected by the global minimum (or 
maximum) of some local driving force. The activity of 
the model in the critical steady state is highly correlated. 
It forms an anisotropic fractal in d + 1 (d spatial and one 
temporal) dimensional space |40|. An example of this 



fractal activity pattern for the one dimensional evolution 
model was shown in Fig. |l|. We recall that one charac- 
teristic exponent of this fractal is its mass dimension D. 
This dimension relates S 1 , the total amount of activity 
within a certain time interval to the spatial extent, or 
range of activity, R, through 



R 



(46) 



Due to the sequential character of activity, S is trivially 
equal to the number of time steps within a selected time 
interval. 
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We have already mentioned that in extremal models 
there exists a purely geometrical fractal property of the 
pattern of activity. The distribution of distances (Levy 
jumps) between subsequent minimal sites, f m in(s), obeys 
a power law Pj um p{r) ~ r _7r . This fractal property was 
described using backward avalanches in Section IVE. 

The pattern of activity has another feature where cuts 
in different directions are fractals themselves. Since the 
pattern is anisotropic, cuts in different directions have 
different mass dimensions. The cuts in the spatial direc- 
tion, at an arbitrary point in time, s, are in fact trivial: 
there is only one active site at any given time. Each 
cut has exactly one point of activity. The fractal prop- 
erties of cuts in the time direction, s, at a given point 
in space, are more interesting. Looking at Fig. El one 
can see that the activity has a tendency to revisit sites. 
At any given site, the activity is recurrent in time, and 
can be considered to be a "fractal renewal process" [(75) . 
The collection of return points on the one dimensional 
time axis forms fractal with dimension < d < 1. As 
in earlier sections, we assume that the projection of the 
activity pattern onto the original d dimensional lattice 
forms a dense, compact region of dimension d. Thus we 
can consider the activity pattern within a spatial region 
R d , to be composed of R d one dimensional fractal time 
lines. In order to encompass all of the activity, each time 
line has a length ~ R D . The number of points of the ac- 
tivity cluster which fall on any given time line scales as 
d. Here, we assume also that there is only one dimension 
d and no multifractal properties of these points on each 
time line. Thus the quantity R d R Dd , the total number of 
active sites in all time lines, must scale in the same way 
as the mass of the entire cluster, S ~ R D . This gives an 
exponent relation: 

d = 1 - cl/D . (47) 

Note that the exponent d is not defined for IP. 

Complex spatio-temporal fractal patterns also can be 
observed in non-SOC systems when their parameters are 
fine tuned so that they become critical. These avalanche 
patterns again are characterized by the fractal dimen- 
sions of different cuts. Exponent relations analogous to 
Eq. (|47]) are obtained below. 

For example, let us consider a large, finite directed per- 
colation cluster on a d + 1 dimensional lattice. A part of 
such a cluster is shown in Fig. ^9| for d = 1. This cluster 
is asymmetric with respect to the t direction. Recall that 
in DP all sites are updated in parallel, so that t — ► t + 1 
when all n(t) sites have been advanced. Self-similarity 
requires that the duration T scales with the spatial ex- 
tent in any one of the d directions perpendicular to time, 
R, as T ~ R z where z is the usual dynamical exponent 
relating space and time. The total size of the cluster, S, 
scales with the spatial extent as S ~ R D , where D is the 
avalanche dimension. Unlike the previous sequential ex- 
ample, we now have more than one active site at a given 



time step. Thus, T and S represent different quantities 
and the corresponding exponents z and D, relating them 
to spatial size R are not equal to each other. In order 
to compute D, usually the cluster is partitioned into R z 
equal time slices. Each such slice contains n act ~ R daat 
points of the cluster. This method of partitioning the 
cluster gives R D ~ R z R dact , so the avalanche dimension 
is given by D = z + d act . As in the previous example, the 
avalanche cluster can be composed as R d one-dimensional 
fractals, parallel to the time axis, which each contain T d 
parts of the cluster. Consequently, 

R D ~ R d R zd ■ d = ^—^ . (48) 

Note that this relation contains Eq. (f47j ) as a limiting 
case when z = D and d act — (one active site at any 
given time step). 

Knowledge of the exponent d enables us to calculate 
two important distributions characterizing return times 
of activity to a given point in space. As a result, it also 
characterizes the power spectrum of local activity. The 
first return probability distribution Pfirst (t) is the dis- 
tribution of 'hole' sizes, or intervals, separating subse- 
quent return points of activity. Here t is the size of the 
hole (either in parallel time t or sequential time s). This 
distribution is normalizable; J Q PFiRST\t)dt = 1. The 
average total number of return points, n(T), in an inter- 
val of length T is given by the fractal dimension of return 
points as n(T) ~T d . It can be related to the first return 
probability; 

n{T) = T — n(T) J P FIRST {t)tdi , where 

P F iRST{t) ~ r TF '» ST for t » 1 . (49) 

If tfirst < 2 then the divergence at the upper limit 
must cancel the T term, so that T ~ n{T)T 2 ~ TFIRST . 
This leads to the scaling relation, 

d = t~first - 1 , (50) 

connecting the fractal dimension of return points to the 
distribution of hole sizes. 

The second distribution, Pall [t-i t) is the probability 
that activity at position at time will be at r at time 
t. This quantity does not obey the same normalization 
condition as Pfirst- Instead J PALL(r,i)df — N, where 
N is the average number of active sites. -Pall(0,?) is 
the probability for the activity at time t to revisit a site 
that was visited at time 0. Unlike Pfirst (i) it doesn't 
require that the return is the first return of activity, so it 
is often referred to as the distribution of all return times. 
Since n(T) is simply the sum of all returns of activity to 
a particular site up to time T, we have 

n(t+ 1) - n(t) = P A LL(0,t) . (51) 
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Assuming a power law form Pall(0, t) ~ t TALL for i 
1, one gets for tall 



d = 1 — TALL 



(52) 



Comparing equations ( |50| ) and (|5^) gives the general re- 
lation 



TFIRST + TALL — 2 



for 



TFIRST < 2 



(53) 



connecting the "lifetime" exponents for the first and all 
returns of activity. 

Recently, Ito J7(| has examined the International Seis- 
mological Center data of California earthquakes and has 
found that the "all" and "first" return time distributions 
for earthquakes to return to a given location are power 
laws over approximately two decades with characteristic 
exponents that obey Eq. (53); i.e. tfirst ~ 1.4 and 
t all ~ 0.5. In addition, he found that the distribution 
of jumps between subsequent hypocenters of earthquakes 
may also be a power law with an exponent tt w 1.7. This 
is a remarkable demonstration of the generality of our 
results. 



A. 1/f Noise and Punctuated Equilibria 

Since Pall(0,£) is the autocorrelation function of the 
activity, the power spectrum is simply 



S(J) 



+ 00 



P A LL(0,t)e 2 " ifi dt ~ i 
J 



(54) 



The mathematical relationship between return times and 
the power spectrum was derived previously using differ- 
ent methods by Lowen and Teich J75| . Here we have 
shown that 1/f type noise emerges naturally in both 
self-organized and non-self-organized critical systems as 
a consequence of avalanche dynamics. As a result, the ac- 
tual exponent d that characterizes the noise is determined 
by the dimension D of the avalanches. Eqs. (^?j), (48), 
([si]) establish a formal connection between 1 // noise and 
fractal scaling behavior, i.e. spatio-temporal complexity, 
in both tuned and self-organized critical models. 

Model dependent behavior occurs between the upper 
critical and lower critical dimension. In the mean field 
limit, or above the upper critical dimension, the activity 
is barely able to return and tfirst = tall = 1. As 
a result, the power spectrum, S(f) ~ 1//°, corresponds 
to white noise. On the other hand, at the lower critical 
dimension, the activity becomes dense in time and d 
1. In this case, the power spectrum S(f) ~ 1//, with 
logarithmic corrections. 

Eqs. (^7] - were checked by numerical simulations. 
We simulated bond DP on a square lattice in 1+1 di- 
mensions at / = 0.6445 for L = 3000. The data shown 
in Fig. [50[ tfirst — 1-86 and tall — 0.14, are in agree- 
ment with the theoretical prediction tfirst = 1.84±0.02 



and tall = 0.16 ± 0.02 based on the exponents D and 
z in Ref. @ . Also S(f) ~ l// 084 in 1+1 dimen- 
sions. In d = 1 we simulated the BS branching process 
at branching probability / = 0.667 and averaged over 
w 10 9 mutations to obtain Fig. [Si]. Our measured values 
are t F irst = 1-58 ± 0.02 and tall = 0.42 ± 0.02, quite 
close to the predicted values 1.59 and 0.41, respectively. 
Similar results were found for d — 2, at branching prob- 
ability / = 0.390, tfirst ^ 1-28 and tall ^ 0.70 to be 
compared with 1.31 and 0.69 from the formulae above. 
The predicted power spectrum is S(f) ~ i/j - 59 m d = 1 
and S(f) ~ l// ' 31 in d = 2. We measured the power 
spectrum S(f) ~ 1/f - 58 for the one dimensional evolu- 
tion model, as shown in Fig. |32]a, and S(f) ~ l// 31 for 
the two dimensional evolution model, as shown in Fig. 



The evolution model was introduced in an attempt to 
model punctuated equilibria in biological evolution. Fig. 
|33] shows the accumulated number of activations, or re- 
turns to a given site, for the one dimensional evolution 
model. In terms of the model, these accumulated re- 
turns would roughly correspond to accumulated muta- 
tions in a given species. The resulting devil's staircase 
shows plateaus of stasis interrupted by bursts of activity. 
The plateaus have a power law distribution in sizes decay- 
ing as S~ TFIRST . This staircase behavior is qualitatively 
similar to the punctuated equilibrium behavior observed 
for the evolution of real species [Q. This suggests that 
the punctuations for a single species are correlated to the 
avalanches in the global ecology. 



VI. INTERFACE DEPINNING 

So far, the main emphasis in this article has been 
on a class of extremal SOC models. In chapter IV, we 
started comparing these extremal models with their con- 
stant force parallel counterparts, where all currently ac- 
tive sites are updated in parallel at each time step. Here 
we discuss the relationship between tuned and SOC mod- 
els of interface depinning in more detail. 

The behavior of an interface driven in the presence 
of quenched random pinning forces appears in a wide 
variety of contexts. These include, among others, fluid 
invasion in a porous medium (59), the motion of mag- 
netic domain walls Jm| , flux lines, or charge density waves 
|]l7f , [7^1 , {ft! in the presence of quenched disorder. The 
depinning transition in the charge density wave system 
has previously been described in terms of avalanche dy- 
namics in sandpile models ]79|] . The notion of an 
external force can be easily incorporated into the rules 
of the models we have discussed. Applying an external 
force f means moving in parallel all sites having their 
local "depinning" force fa < f . Below the threshold 
force f c , the interface is pinned in one of many possible 
metastable states. An external force greater than the 
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threshold causes the interface to move with a finite aver- 
age velocity which vanishes continuously at a depinning 
transition. The instantaneous velocity proportional to 
the density of active sites, fluctuates in time around its 
average value. 

Instead of tuning the force to the depinning thresh- 
old, the velocity may be tuned. In this case, the force 
fluctuates in time, in such a way that the instantaneous 
velocity is constant. The depinning transition at con- 
stant velocity is reached as the limit when the imposed 
velocity vanishes. The instantaneous interface velocity 
in discrete lattice models is simply proportional to the 
current density of active sites: v = cn ac t/L d , where c 
is a constant of 0(1). We restore the rules of SOC ex- 
tremal interface models if we require that strictly one site 
is moving at any given time step and this site is selected 
as most unstable in the whole system, i.e the global min- 
imum of /j. This corresponds to the interface velocity 
v e xt = cL~ d . Some of the critical properties of the model 
driven at constant velocity are the same as at constant 
force but others are different. It is interesting to note that 
Wilkensen and Willemsen ]33| introduced invasion per- 
colation as a modification of earlier models |}4| that were 
driven at constant force. As they noted, the introduction 
of extremal dynamics corresponds to the constant veloc- 
ity invasion process in the limit of vanishing velocity. 

Theoretical studies of interface depinning have concen- 
trated along two main directions. One theoretical ap- 
proach has been to apply a functional renormalization 
group procedure to a variety of different depinning phe- 
nomena j7j|, pC| . This has yielded both perturbative 
results and results that have been claimed to be exact. 
For example, for the LIM defined by Eq. (jl|), it has been 
claimed 65 that the roughness exponent x = (4 — d)/3 



exactly. Another theoretical approach has been to con- 
sider simple lattice models for these phenomena, which 
may, due to their simplicity, be analytically tractable. 

Several of the parallel lattice models defined in Section 
II describe interface depinning at constant force. The 
LIM and TLB models presumably describe the depin- 
ning of magnetic domain walls having purely linear force 
terms, and other elastic interfaces having nonlinear force 
terms, respectively [ jBl| . Sneppen was the first to intro- 
duce a self-organized critical extremal model dedicated 
to interface depinning; similarly an SOC, or constant ve- 
locity, variant of the LIM model can be constructed jlO| . 
The general properties of these SOC models have been 
analyzed in the preceding sections. Our main results were 
the derivation of two exact equations, a stationarity con- 
dition, and numerous exponent relations. These scaling 
relations are summarized in Table I. 

We now discuss how our results may be generalized to 
the tuned depinning models driven by constant external 
force. Here, we only analyze the behavior as the depin- 
ning transition is approached from below. In this limit, 
we will show that the tuned and SOC variants have the 



same avalanche dimension D, the same avalanche dis- 
tribution exponent r, and the same roughness exponent 
X- However, the exponents that describe propagating ac- 
tivity within avalanches are different. In particular, the 
stationarity law r\ — does not necessarily hold for the 
tuned models at constant force, although it does hold for 
the extremal SOC models. This implies that the fractal 
dimension of activity df and the dynamical exponent z 
can also be different in these two cases. These results are 
verified numerically. 

We first concentrate on the case of the Sneppen and 
TLB models which are constant velocity and constant 
force versions of the same depinning phenomenon. Let us 
recall the definition of the f a - avalanche in the Sneppen 
model. This avalanche intervenes between subsequent 
punctuations of the barrier f by the signal / m in(s), 
which is the extremal value of the pinning force. Geo- 
metrically, an f a - avalanche takes the interface from one 
critical "blocking surface" where all the random numbers 
are greater than or equal to f to another f a blocking 
surface, as shown in Fig. ||. Tang and Leschhorn p4| ] 
showed that in one dimension, these blocking surfaces 
correspond to percolating paths on a directed percola- 
tion cluster, formed by all sites with < f Q . Given a 
collection of random pinning forces f(x, h), and an initial 
interface configuration identifying with a critical blocking 
surface, the next blocking surface that is encountered un- 
der Sneppen dynamics is the same blocking surface that 
would be encountered in an equivalent parallel model, i.e. 
the TLB model, driven with an applied force f . The 
order in which the sites between these two blocking sur- 
faces are invaded is completely different, though. In the 
Sneppen model, one always chooses the smallest random 
pinning force, f m in- In the TLB model, one advances all 
sites where fi < f in parallel. Nevertheless, the differ- 
ence between the initial interface configuration and the 
final configuration, given by the two blocking surfaces, is 
the same for the two models. Thus the TLB model has 
the same threshold f c , roughness exponent x of blocking 
surfaces, and avalanche distribution P(S,f ). This has 
been confirmed numerically in Refs. p2], [Q, |32|, fl45| , 

& 

The roughness exponent x characterizes the saturation 
width, w, of the height fluctuations in a system of size 
L, so that w 2 =< (h- < h >) 2 >~ L 2x f|, §§. The 
roughness exponent x an d the avalanche dimension D 
are related via 



D = d + X 



(55) 



The logic behind this relation is: 1) The set of sites ad- 
vanced at least once during the course of an avalanche 
is assumed to form a compact object having the same 
dimension d as the interface substrate. In d = 1 any 
connected set is also compact so this assumption is ob- 
viously valid. In higher dimensions, we are not aware 
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of any proof that the avalanches are compact. 2) It is 
assumed that there are no multifractal features of the 
measure, defined as a total number of updates during an 
avalanche, on the set of n cov sites. Namely, there is only 
one scale characterizing the number of times each site 
covered by the avalanche is updated. If both conditions 
1) and 2) are satisfied, the avalanche volume, which in 
extremal models is simply proportional to its temporal 
duration S, can be written as R d R x , and the relation 
( p5| ) is satisfied. On the other hand, it is not completely 
inconceivable that in high enough dimensions one of the 
assumptions can be wrong for the interface models. In 
particular, this would occur in dimensions above the up- 
per critical dimension of the model, if an upper critical 
dimension exists. Then Eq. (|55| ) no longer holds, and 
D, the avalanche dimension, can be smaller than d, the 
dimension of the substrate. The results of the numerical 
simulations of the TLB model reported in |44j , 0] sug- 
gest that the relation ( |55| ) may be valid in dimensions as 
high as d = 4. 

Taking into account these assumptions, our results 
contained in section IV B, which were derived for Snep- 
pen model directly apply to its parallel counterpart - the 
TLB model. 

The critical dynamics of the interface as it propagates 
from one blocking surface to another may be different 
in the constant velocity (SOC) and constant force inter- 
face models. This difference is not restricted to a trivial 
time redefinition, where in extremal versions the tempo- 
ral duration of an avalanche is simply proportional to its 
volume, while in parallel versions these quantities differ 
because all active sites are advanced in parallel. What is 
more interesting is that even the fractal dimension, d/, 
of active sites (sites with fa < f ) is different in these 
two cases. Since the dynamics of the Sneppen model is 
stationary, the dynamical exponent t] = 0. There is no 
such stationarity condition in the TLB model at constant 
force; thus r) is not necessarily zero in that model. In fact, 
in one dimension Tang and Leschhorn |62| have argued 
that z = 1 exactly for TLB at constant force, in contrast 
to z = D — df ~ 1.21 which is predicted as a corollary 
of the 77 = law for the Sneppen model and which was 
measured in Fig. In fact, we measured z for the 

TLB model, as also shown in Fig. 
z = 1 prediction and previous numerical results 



and confirmed the 

HID- 

This demonstrates that the dynamical critical exponents 
in SOC and tuned critical phenomena can be different. 

It is possible, though, that these differences may be- 
come smaller or disappear entirely for tuned and SOC 
versions as the substrate dimension increases. The nu- 
merical simulations of the TLB model in dimensions up 
to d = 4 H], Ef| may possibly support this point of 
view. One of the consequences of the 77 = law is that 
the exponents they defined as r surv and 8 would sat- 
isfy the relation 8 = df/z — r surv — 1. This relation is 
obviously violated in the one dimensional model where 



they @ report r£+ 1} = 1.46(2) and = 0.60(3), 

which once again rules out 77 = for the one dimen- 
sional TLB model. Their results for higher dimensions 
<5< 2+1 ) = 1.14(6), t£XV = 2.18(3), <5( 3+1 > = 1.6(1), 



.(3+1) 



2.54(5), <5< 4+1 ) = 1.9(2), 



.(4+1) _ 



3.0(2), how- 



ever, seem to obey the r\ = relations within their numer- 
ical uncertainties. We do not know whether this appar- 
ent agreement in dimension d > 1 is simply a numerical 
coincidence. 

A similar argument that in tuned and self-organized 
versions of the model the exponents r and D are the 
same, while rj in general is different can be applied to 
the LIM model. In this model, the advancement at any 
given site can never cause a neighboring unstable site to 
become stable; the motion at an active site can never 
destroy another active site. This means that one can 
interchange, arbitrarily, the order in which the unstable 
sites move without changing the final metastable config- 
uration that is reached. This property is reminiscent of 
the Abelian properties of the sandpile model exploited by 
Dhar |2j| . This interchangeability means that the critical 
force f c , as well as the exponents D and r describing the 
avalanche statistics, are the same for the two versions. 
Obviously, though, r\ can be different because the dy- 
namics during an avalanche depends on the way unstable 
sites move. In fact, Leschhorn measured z = 1.42 ± 0.03 
in one dimension and z = 1.58 ± 0.04 in two dimension 
for the tuned LIM at constant force; whereas, based on 
our measured values of D and r, we predict quite differ- 
ent values z ~ 1.94 (one dimension) and z ~ 2.04 (two 
dimensions) for the SOC or constant velocity variant. 

For the interface models, it is conventional to assume 
the validity of the relation (|55|) and to express the critical 
exponents in terms of \ — D — d, and is, the spatial 
correlation length exponent along the substrate. In this 
case, Eqs. ( |3"Tf - ^2| ) can be rewritten in terms of v and x 
as 



r = 1 



d- l/v 



X 



1 = 1 + XV 



(56) 



(57) 



In the one-dimensional Sneppen model, v and \ 
have been derived from the fractal properties of the 
1+1-dimensional directed percolation cluster pil, p5[ . 
Namely x = v® 1 * JvP p and v = i>P p , where vx p and 



r 



11 



are parallel and perpendicular correlation length ex- 
ponents in 1+1-dimensional directed percolation. For the 
one dimensional Sneppen model, these expressions dif- 
fer from the predictions of Olami, Procaccia, and Zaitek 
that t = 2/(2 + x) and 7 = 2 g. The avalanche di- 
mension D was measured for the self-organized LIM as 
shown in Fig. |l|. In d = 1, D = 2.23 ± 0.03, and in 
d = 2, D = 2.725 ± 0.02. For the LIM in both one 
and two dimensions, D was measured by computing, for 
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each backward avalanche, the duration S of the back- 
ward avalanche, and the squared distance R 2 from the 
site starting the backward avalanche to the site ending 
it. Substituting these measured values into the relation 
X = D — d, gives x — 1-23 in d = 1 and x — 0.72 in 
d = 2. These values for x are i n agreement with numer- 
ical simulations by Leschhorn [ p8| , but higher than the 
prediction x — (4 — d)/3 from functional rcnormaliza- 
tion group (RG) calculations |35| in both one and two 
dimensions. 

One might wonder if the one dimensional LIM de- 
scribes the behavior of any real physical system. Our 
measured value of the interface roughness x — 1.23 ±0.02 
derived from the avalanche fractal dimension D = 
1 + X, as well as the reported values x — 1-25 ± 0.01 J||] 
and x — 1.2 ± 0.1 ^6|, contradicts the usual condition 
that the " self-affine" interface looks flat when viewed at 
a sufficiently large length scale. It may be more appro- 
priate to call the interfaces with x > 1 "super-rough" 
|H , Q . As was suggested in Q , one possible physical 
realization of the super-rough interface is fluid invasion 
in a porous medium. Martys, Robbins, and Cieplak ]6C| ] 
introduced an explicit model for this process. They have 
shown that depending on the wetting properties of the in- 
vading fluid in 1+1 dimensions one observes two distinct 
universality classes of the constant pressure, i.e. force, 
interface depinning transition. One was identified with 
the constant pressure version of invasion percolation |}4| , 
having an extremely curved interface with overhangs at 
all length scales. The interface from the other univer- 
sality class has overhangs only at small length scales, 
and Roux and Hansen have argued that the 1+1 LIM 
shares the same universality class [ p6[ . Our numeri- 
cal results support this conjecture. We have measured 
r = 1.13 ± 0.03 in very good agreement with the mea- 
sured value r = 1.125 ± 0.025 of (6Q|. Scaling relations 
analogous to Eqs. (7, 31, 32) for the non SOC model 
were derived in [ |60f based on the assumption that the 
avalanches in the model are isotropic; specifically they 
assumed D = d + 1, which is different than our result 
D = d + x- Nevertheless, these authors considered the 
interface itself to be self-affine (x 7^ 1), and not isotropic. 
Within numerical uncertainty, it appears to us that their 
results are also consistent with anisotropic super-rough 
avalanches having the same fractal properties as the in- 
terface itself (x = D — d). For example, inserting their 
numerical value v = 1.30 ± 0.05 and our numerical value 
D = 2.23 into Eq. ( p36| ) one gets t = 1.105, which is also 
consistent with the measured values cited above. The 
question of the shape of the avalanches in their model 
could probably be resolved by direct measurements of 
the avalanche volume S vs. its spatial extent R. 

Similiar considerations may apply to self-affine inter- 
faces in the 2+1-dimensional Random Field Ising Model 
studied by Ji and Robbins plf . They measured r = 
1.28 ± 0.05 and v = 0.75 ± 0.05 in agreement with our 



numerical result r = 1.29 + 0.02 for the 2+1-dimensional 
LIM and Eq. @ using D = 2 + x = 2.72 + 0.02. Again, 
they assumed the avalanches to be isotropic in the Ran- 
dom Field Ising Model, while the avalanches in the LIM 
are anisotropic and self-affine. 



VII. CONCLUSIONS 

We have presented a comprehensive theory for 
avalanche dynamics in evolution, growth, and depinning 
models. These models are defined in Section II and rep- 
resent different universality classes. We have shown that 
avalanche dynamics leads to spatio-temporal complexity, 
and emerges as a result of extremal dynamics in driven 
systems. Spatio-temporal complexity is manifested in 
the formation of fractal structures, the appearance of 
1// type noise, diffusion with anomalous Hurst expo- 
nents, Levy flights, and punctuated equilibrium behav- 
ior. These phenomena can all be related to the same 
underlying avalanche dynamics. 

We present two exact equations for these phenomena. 
(1) The approach to the critical attractor is governed by a 
"gap" equation for the divergence of avalanche sizes. (2) 
The average number of sites covered by the avalanche 
can be related to the average size of avalanches. If there 
is a power law divergence for the average avalanche size, 
then the number of sites covered by an avalanche diverges 
with exponent -1. In addition, the conservation of activ- 
ity in the stationary state manifests itself through the 
fundamental relation r/ = 0. It follows that many of the 
critical exponents in a class of SOC extremal models can 
be derived from two basic exponents. 

These exponent relations are summarized in Table I. 
Depending on the model it may be more convenient to 
use one or another basic set of two exponents. In Table I 
we write our exponent relations in terms of two such basic 
sets: (D, v), or (D, r). The scaling relations are defined 
for the Bak-Sneppen evolution model, the Sneppen model 
for SOC interface depinning, the Zaitsev model flux creep 
model, an SOC "linear" interface model, and invasion 
percolation. The horizontal lines separate results from 
different sections of this article. 

Our results from numerical simulations are summa- 
rized in Table II. This Table contains the expressions 
for the exponents of the one and two dimensional Bak- 
Sneppen model, one and two dimensional LIM, and the 
one dimensional Sneppen model based on direct numeri- 
cal simulations. The overall consistency between the ex- 
ponent relations from Table I with the numerical results 
in Table II demonstrates the validity of our new scaling 
relations. Most of the numerical results from Table II 
were illustrated in figures throughout the text. 

Some of these new scaling relations show that l/f 
noise and spatial fractal behavior can be unified, and 
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have a natural explanation in terms of avalanche dynam- 
ics in both SOC and non-SOC systems. The pattern of 
avalanches in the critical state can be described as a frac- 
tal in d spatial plus one temporal dimension, with mass 
dimension D. Temporal behavior, such as 1// noise in 
the power spectrum, and spatial long-range correlations 
can be formally related as different cuts in the same 
underlying fractal. In the stationary state, time rever- 
sal symmetry is broken, so that forward and backward 
avalanches have different statistics. We derive a scaling 
relation for the Levy distribution of jumps in the SOC 
extremal models models. Finally, we point out the sim- 
ilarities and differences between interface depinning at 
constant velocity (SOC) and constant force. 
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VIII. APPENDIX A: INVASION PERCOLATION 

Invasion percolation is very sensitive to the definition 
of the boundary of the invaded cluster. In one version 
of the model, the boundary consists of any site having 
at least one nearest neighbor in the invaded cluster. The 
boundary defined by this rule includes sites on the exter- 
nal perimeter of the invaded cluster as well as sites on 
the perimeter of numerous " lakes" in the interior of the 
cluster. The invaded cluster in this model from time to 
time exactly identifies itself with the infinite cluster of 
ordinary percolation. At these instances, all of the ran- 
dom numbers on the boundary are uniformly distributed 
above the "gap" , f c , where f c is the critical density of or- 
dinary percolation. In two dimensions, the whole bound- 
ary of the cluster of ordinary percolation has a fractal 
dimension, d B , equal to the fractal dimension D of the 
cluster itself, and D = d B = 91/48 ~ 1.89 @. We will 
refer to the invasion percolation model with the bound- 
ary defined in this manner as IP-1. 

An important physical realization of invasion perco- 
lation is the displacement of one fluid by another in a 
porous medium. In this case, once a region of non- 
invaded fluid is completely surrounded by the invading 
fluid, no further invasion can take place on this part of the 
boundary due to incompressibility. This event, known as 
self-trapping, can be taken into account in the invasion 
percolation model by changing the definition of the active 
boundary. With self-trapping, only the sites touching the 
external perimeter of the invaded region comprise the ac- 
tive boundary. It is believed that self-trapping does not 
change the fractal dimension of the invaded cluster in 



two dimensions |l4[] . In ordinary two dimensional perco- 
lation, it is known that the external perimeter of the 
cluster has a smaller fractal dimension than the clus- 
ter itself, i.e. d B < D. Depending on the details of 
the precise definition of the external perimeter, one gets 
d B = 7/4 = 1.75 [g, or d B = 4/3 ~ 1.33 ||. These 
different definitions thus give different variants of the in- 
vasion percolation model, which we will refer to as IP-2, 
and IP-3, respectively. 

In Section III we used the result by Roux and Guyon 
po| for the scaling of the boundary of the growing IP 
cluster with its volume. For the sake of completeness, 
we reproduce their arguments. During the transient, the 
invaded cluster can be characterized by a growing cor- 
relation length, £. Since the growth starts from a d- 
dimensional base of the d+ 1-dimensional hypercube, the 
natural scaling form for the invaded volume s is 

s ~ {L/O d e , (58) 

and for the boundary of invaded cluster, b(s), 

b(s)^(L/O d i dB ■ (59) 

Combining these two equations we get: 

b(s)/L d = ( S /L d )^ , (60) 

and therefore g = f,Zt - 

It is interesting to note that the simplest assumption 
that the effective distance from the critical point scales 
as f c - G(s) ~ f-V" ~ {s/L^- 1 /^ "^ agrees with the 
gap equation after we substitute into it the expression 
for g and for 7. This once more confirms the overall 
consistency of our scaling relations. 
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FIG. 1. The space-time fractal activity pattern for the one 
dimensional evolution model. Time is measured as the num- 
ber of update steps, S. 

FIG. 2. Snapshot of the stationary state for the evolution 
model in a finite one dimensional system. Except for a lo- 
calized region, the avalanche, where there are small random 
numbers, all the random numbers in the system have values 
above the self-organized threshold / c = 0.66702. 

FIG. 3. Schematic picture of an avalanche separating two 
interface configurations in the Sneppen model. The size, S, 
of the avalanche corresponds to the shaded area. 



FIG. 5. Illustration of the hierarchical structure of the / - 
avalanches. The big avalanche is subdivided into smaller ones 
as the auxiliary parameter f is lowered. The lines span Jo- 
avalanches for f — 0.63, 0.59, and 0.54, respectively. 

FIG. 6. Distribution of avalanches from simulation of the 
BS branching process in a) one and b) two dimensions. The 
asymptotic slope of the log- log plot gives r = 1.07 ± 0.01 in 
Id and r = 1.245 ± 0.01 in 2d. 



FIG. 7. One realization of the extremal signal fmin as a 
function of s in the stationary state. If the auxiliary param- 
eter f is raised by an infinitesimally small amount df , the 
breaking points that had f m in(s) between f and f ± df 
(filled circles) will no longer stop the f ± df - avalanches, 
and the average avalanche size will increase. 

FIG. 8. The gamma plot (1 — f a )/ < n cov > vs. f for 
the two dimensional evolution model. The intersection with 
the horizontal axis gives / c , and the slope gives I/7. We find 
f c = 0.328855 ± 0.000004 and 7 = 1.70 ± 0.01. 



FIG. 9. As Fig. 8 for the Id evolution model. We find 
f c = 0.66702 ± 0.00003 and 7 = 2.70 ± 0.01. 



FIG. 10. The plot of (1 — f )/ < n cov >/ o as a function 
of f for the one dimensional Sneppen model. The slope of 
the straight line is I/7 = 0.470 and the crossing point with 
/ -axis correspond to / c = 0.465. 



FIG. 11. The average size of avalanches, < S >, vs. (/ c — /) 
for the 2d evolution model. The asymptotic slopes yields 
7 = 1.69 ±0.03. 



FIG. 12. The average size of avalanches, < 5* >, vs. (f c — f) 
for the Id evolution model. The asymptotic slopes yields 
7 = 2.71 ±0.02. 

FIG. 13. a) R vs. S for the 2d evolution model (av- 
eraged over 107000 avalanches). b) Derivative of this 
curve. The asymptotic value corresponds to 1/D and gives 
D = 2.92 ±0.02. 



FIG. 14. R vs. S for the Id evolution model involving 
3 x 10 11 updates. From this data we extract D = 2.43 ±0.01. 



FIG. 4. Value of the extremal signal f m i„ as a function 
of s during the transient in a small (L=20) one dimensional 
evolution model. The full line shows the gap, G(s), as a 
step-wise increasing function of s. The gap is an envelope 
function that tracks the peaks in f m in- 
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FIG. 15. S vs. R from backward avalanches in Id (open 
circles) and 2d (filled circles) LIM. In Id, the slope is 
D = 2.23 ± 0.03 from simulations of 5 x 10 7 time steps in 
the stationary state of a system of size L — 3000. In 2d, the 
slope is D = 2.725 ± 0.02 from simulations of 10 9 time steps 
in the stationary state of a system of size 300 x 300. 

FIG. 16. Typical n(s) for the two dimensional evolution 
model. 

FIG. 17. Average number of active sites < n(S) > vs S 
(lower curve) and average activity of surviving avalanches, 
nsuRv(S) (upper curve), for the two dimensional evolution 
model at f c = 0.3289. The slopes yield d s = 0.25 ± 0.005 and 
t) = 0.0 ± 0.002 respectively. 

FIG. 18. Average number of active sites < n(S) > and 
average activity of surviving avalanches, usurv(s), for the 
one dimensional evolution model at f c = 0.6670. The slopes 
yield d s = 0.11 ± 0.02 and 77 = 0.01 ± 0.02. 

FIG. 19. Parallel time, t, vs. spatial extension, r, for 
the one dimensional Sneppen model (open dots) and the 
one dimensional TLB model (filled dots). The slopes gives 
2 = 1.19 ± 0.05 for the Sneppen model and 2 = 1.00 ± 0.05 
for the TLB model. 

FIG. 20. Parallel time, t, vs. spatial extension, r, for a 
the two dimensional (small filled dots) and one dimensional 
evolution (large open dots) models. The asymptotic slopes 
are 2 = 2.17 ± 0.02 and 2 = 2.10 ± 0.05, respectively. 

FIG. 21. < n(S) > and nsunv(S) for the one dimensional 
Zaitsev model with L = 10 4 . d s = 0.15 ± 0.03. The asymp- 
totic value of the slope of < n(S) > is r\ — 0.0 ± 0.006. 

FIG. 22. The hierarchical structure of backward and for- 
ward avalanches. The avalanches that started inside a larger 
parental avalanche are completely contained within it. 

FIG. 23. The overall distribution of backward avalanches 
in the Id Sneppen model. The slope of the curve gives 
r 6 a " = 3 - t = 1.745 ± 0.02. We simulated 10 7 time steps 
in a system of size L = 1000. 

FIG. 24. The overall distribution of backward avalanches 
in Id (open circles) and 2d (closed circles) LIM. In Id, the 
slope of the curve corresponds to r b aii = 3 — r = 1.87 from 
simulations of 5 x 10 7 time steps in a system of size L — 3000. 
In 2d, the slope of the curve corresponds to t£ u — 3 — r = 1.71 
from simulations of 10 9 time steps in a system of size 300 x 300. 



FIG. 25. Distribution of avalanches from simulations of the 
Id Sneppen model for L = 2 x 10 5 and f c = 0.4614. The 
asymptotic slope of the log-log plot gives r = 1.26 ± 0.01. 

FIG. 26. The distribution of jumps between subsequent 
minimal sites in the Id evolution model. The slope of the 
straight line corresponds to tt — 3.23. We simulated 5 x 10 7 
time steps in the system of size L = 3000. 

FIG. 27. The distribution of jumps between subsequent 
minimal sites in the Id LIM. The slope of the straight line 
is ty = 2.93. We simulated 5 x 10 7 time steps in the system of 
size L = 3000. 

FIG. 28. The distribution of jumps between subsequent 
minimal sites in the 2d LIM. The slope of the straight line 
corresponds to ir — 2.89. We simulated 10 9 time steps in a 
system of size 300 x 300. 

FIG. 29. The pattern of active sites in 1+1-dimensional 
bond directed percolation. 

FIG. 30. Distribution of return times for bond DP on a 
square lattice for / = 0.6445 and L = 3000. The asymptotic 
slopes give tfirst — 1-86 and tall — 0.14. 

FIG. 31. Return times for the a) Id BS branching process 
at branching probability / = 0.667 averaged over fa 10 9 mu- 
tations. The corresponding exponents are tfirst — 1-58 and 
tall — 0.42. b) 2d BS branching process at / = 0.3289 aver- 
aged over « 10 9 mutations, tfirst — 1-28 and tall — 0.70. 

FIG. 32. Power spectrum for a) the one and b) the two 
dimensional evolution model. 

FIG. 33. Punctuated equilibria. The curve shows the accu- 
mulated number of times a specific site is visited, or is the min- 
imal site, in the one dimensional evolution model. This would 
roughly correspond to the accumulated number of mutations 
in a particular species. The pattern of change is step-wise 
rather than being smooth. 
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TABLE I. 


Exponent relations for the extremal models in d dimensions. 


For invasion percolation, d should be replaced with 


the fractal dimension of the active boundary ds except in the scaling relation for p (see text). 


The exponents tall, tfirst, 


and d are not defined for invasion percolation. In the interface models D is 


related to the interface roughness \ via D — d + \- 


Exponent 


Physical Property 


(y and D) 


(r and D) 


D 


Avalanche Dimension 


basic exponent 


basic exponent 


V 


Correlation Length 


basic exponent 


l 

d-D(r-l) 


T 


Avalanche Distribution 


1 + d -^ /v 


basic exponent 


1 


Average Avalanche Size 


1 + u(D - d) 


2-T 

±-\-d/ D t 


a 


Avalanche Cutoff 


1/vD 


i + i-T 


n 


Relaxation to Attractor 


i 

v(D-d) 


1-d/D 


V 


Average Growth of Activity 


u 


U 


(if 

J 


Dimension of Active Sites 


d- i 


D(t - 1) 


Z 


Dynamical Exponent 


D-d+l 


D(2-t) 


IT 


Jump of Minimal Site 


l+D-d+± 


1 + D(2-t) 


all 

T f 


All Forward Avalanches 


2 


2 


all 


All Backward Avalanches 


\-l 


3-r 


TFIRST 


Punctuated Equilibrium 


Z D 


r\ d 

Z D 


TALL 


All Returns 


d 
D 


d 
D 


d 


1// Noise 


1 - ± 
1 D 


1 - ± 
x D 



TABLE II. Critical exponents measured and illustrated in figures throughout the text. All values were determined indepen- 
dently. (• • •) indicates uncertainty in the last digit. Within these uncertainties, all exponents are consistent with the exponent 
relations from Table I. 



Exponent 


Bak-Sneppen model 


Linear Interface model 


Sneppen model 




1 dimension 


2 dimensions 


1 dimension 


2 dimensions 


1 dimension 




f c = 0.66702(3) 


f c = 0.328855(4) 






f c = 0.4614(4) 


D 


2.43(1) 


2.92(2) 


2.23(3) 


2.725(20) 




T 


1.07(1) 


1.245(10) 


1.13(2) 


1.29(2) 


1.26(1) 


7 


2.70(1) 


1.70(1) 






2.13(20) 


71 


3.23(2) 




2.93(3) 
0.000(6) b 
0.15(3) b 


2.89(3) 




n 


0.01(2) 


0.000(2) 




0.03(5) 




0.11(2) 


0.250(5) 




0.26(2) 


z 


2.10(5) 


2.17(2) 






1.19(5) c 


TFIRST 


1.58(2) 


1.28(3) 








TALL 


0.42(2) 
0.58(3) d 


0.70(3) 
0.31(3) d 








d 









b Measured for Id Zaitsev model, presumably from the same universality class |5q ]. 

c This exponent for the TLB model (parallel analog of Sneppen model) was measured to be 1.00 ± 0.05. 
d This exponent was measured from the power spectrum of activity. 
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